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ABSTRACT 


Description, Analysis and Prediction of Sea-—Floor Roughness 


Using Spectral Models 


CHRISTOPHER GENE FOX 


A method has been developed which allows a valid statistical model 
of the variability of oceanic depths to be derived from digital bathy- 
metric soundings. Existing bathymetric contour charts represent low- 
frequency, deterministic models of the sea floor. To describe the 
higher frequency variability, or roughness, of the sea floor requires 
the development of a valid stochastic model. The statistical model of 
sea-floor roughness is also a valuable tool for predicting acoustic 
scattering and in addition contains a wealth of geological information 
for interpreting deep-sea processes. 

To allow the variability of depths to be described as a function of 
scale (spatial frequency), the amplitude spectrum is employed as the 
fundamental statistic underlying the model. Since the validity of the 
amplitude spectrum depends upon the assumption of a statistically sta- 
tionary sample space, a computer algorithm operating in the spatial 
domain was developed which delineates geographic provinces of limited 


statistical heterogeneity. Within these provinces, the spectral model 


is derived by fitting the amplitude estimates with one or two power law 
functions. 

The distribution of the model parameters is examined for a test area 
adjacent to the coast of Oregon. The distribution of roughness corre- 
sponds generally with the various physiographic provinces observed in 
the region. Additional complexities are apparent in the roughness model 
which can not be inferred by simply studying the bathymetry. These pat- 
terns are related to geological processes operating in the region. 

In many cases, the roughness statistics are not constant when 
observed in different directions, due to the anisotropic nature of the 
sea-floor relief. A simple model is developed which describes the 
roughness statistics as a function of azimuth. The parameters of this 
model quantify the anisotropy of the sea floor, allowing insight into 
the directionality of the corresponding relief-forming processes. 
Finally, the model is used to successfully predict the roughness of a 
surface at scales smaller than those resolvable by surface sonar 
systems. The model regression line (derived from a hull-mounted sonar) 
is compared to data from deep-towed sonars and bottom photographs. The 
amplitude of roughness is predicted to within half an order of magnitude 


over five decades of spatial frequency. 
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1. Introduction 


The contour map of oceanic depths, or bathymetric chart, has pro- 
vided a fundamental tool for inferring deep-sea processes, both sedi- 
mentological and tectonic. The identification of the major features of 
the sea floor led to more elaborate geophysical studies and such unify- 
ing theories as sea-floor spreading. The efforts of numerous institu- 
tion. world-wide have produced comprehensive contour charts of global 
bathymetry, which have formed the basis of further geophysical survey 
efforts. 

The method of contouring has the mathematical equivalent of fitting 
a continuous surface to discrete three-dimensional data. In the case of 
bathymetric contours, the discrete data are usually in the form of vari- 
ably spaced and randomly oriented ship tracks along which are discrete 
soundings at some interval. Modern soundings are normally estimated 
from surface ships using acoustic sounders in which the two-way travel 
time of a pulse of sound is interpreted as a depth at a discrete point. 
Geometric spreading of the sound with depth causes a large area of the 
sea floor to be insonified (often called the “footprint”), and the meas- 
ure of the first significant return of sound to the ship is interpreted 
as a sample of the shallowest depth from this large insonified area. 
The uncertainties in the precise location of the measured depth, the 
sound velocity structure of the water column and the positioning of the 
ship at sea, combine to introduce a noise component into the output sig- 
nal. The random noise of measurement, in addition to the uncertainty of 


interpolation between often widely spaced samples and ship tracks, 


imposes a physical limit on both the vertical (depth) accuracy and the 
spatial resolution possible in utilizing bathymetric contouring methods. 

Despite its inherent uncertainties, a bathymetric contour chart 
does provide an absolute estimate of depth at all points in space. In 
the terminology of numerical modelling, a contour chart can be consid- 
ered a “deterministic model” of sea-floor topography. Because such a 
model is analogous to fitting a least-squares surface through randomly 
spaced, noisy data, it is by nature a smoothed representation of the 
actual “submarine topography. The degree of smoothing required (or 
equivalently the cutoff frequency for low-pass filtering of spatial 
frequency) depends upon the accuracy, resolution, and density of data 
used in the contouring. Van Wyckhouse (1973) demonstrated the deter- 
ministic aspect of bathymetric charts with the creation of SYNBAPS 
(Synthetic Bathymetric Profiling System), a data base containing depths 
on an evenly spaced grid. The grid cell spacing of 5° of latitude and 
longitude used for SYNBAPS, seems to represent a workable estimate of a 
suitable interpolation interval for deterministic modelling. Recent 
work at the U.S. Naval Oceanographic Office has extended this 5' grid 
world-wide. 

There are also some practical considerations, in addition to the 
physical limitations, in determining submarine topography to high spa- 
tial frequencies. Survey instruments such as deep-towed sleds incorpo- 
rating stereo photography, narrow-beam sounders and side-scan sonar, are 
able to map small areas of the sea floor with high resolution. However, 
it is practically impossible to extend these surveys globally, due to 
the operational difficulties and the great expense of these methods. 


Another practical difficulty of extending deterministic, gridded models 


to a smaller grid cell spacing is the storage requirements of the data. 
For example, the 5° gridded data set extended from 70°S to 70°N being 
developed by the Naval Oceanographic Office will require approximately 
five million storage locations. To extend this grid to the order of 100 
meters spacing (assuming this was determinable globally) would require 
approximately 5 x 1010 storage locations. 

In light of the physical and practical limits of describing high 
spatial frequency sea-floor topography deterministically, the apparent 
alternative is some probabilistic (stochastic) model describing the dis- 
tribution of features in a statistical sense rather than determining the 
exact locations of depths. To be useful, this stochastic model must 
describe, with a reasonable number of parameters, most of the true vari- 
ability of a region of sea floor. Using a stochastic representation of 
high frequency feature distribution in combination with the lower fre- 
quency deterministic models being developed, an essentially complete 
description of the sea floor is possible. The development of such a 
stochastic model of sea-floor topography is the intent of this 


dissertation. 
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2. Applications 


In addition to the intellectual satisfaction of completing the des- 
eription of sea-floor depths, a stochastic model of high spatial fre- 
quency submarine topography, or surface roughness, has many practical 
scientific and engineering applications. Many of the most direct appli- 
cations are in the field of underwater acoustics. Clay and Medwin 
(1977) provide one of the best physical descriptions of the interaction 
of an acoustic signal with a rough surface. Matthews (1980) reviews the 
importance of relative scale in acoustic bottom interaction. Briefly, 
the relative spatial frequencies of the incident acoustic signal and the 
bottom roughness determine whether the surface acts as a reflector or as 
a scatterer of energy. This relationship illustrates the necessity of 
describing bottom roughness in terms of spatial frequency. 

The importance of acoustics to marine geophysical surveying systems 
cannot be overstated. Sea-floor bottom loss of sonar and seismic sys- 
tems, ranging of side-scan sonars, and signal strength of outer beams on 
multibeam sounders all depend heavily on bottom roughness. In naviga- 
tion applications, the roughness of the sea floor has an impact on the 
backscattering of Doppler sonar, as well as determining the background 
noise for navigation by bottom features. A major application for the 
U.S. Navy of bottom roughness information is as environmental input into 
long-range acoustic propagation models, used in submarine and surface 
ship tracking. For engineering applications, the a priori knowledge of 


the spatial frequency composition of the bottom could be used in 


computer-aided design of future sensing systems, as well as aid in main- 
taining coherent signals for underwater communication systems. 

Another field with significant applications for sea-floor roughness 
information is physical oceanography. This requirement led to the 
extensive work of T. H. Bell. Planetary Rossby waves are strongly 
affected by bottom roughness in long wavelengths (Rhines and Bretherton, 
1973). The propagation of long surface waves such as tides and tsunamis 
Bee also affected (Rhines, 1977). Bell (1973, 1975a) showed that the 
interaction of deep ocean currents and bottom topography may lead to the 
generation of internal gravity waves in the oceans, a major influence in 
ocean dynamics as well as submarine operations. 

Another geophysical application is in the general field of survey 
design. Davis (1974) has formulated a method which, with a knowledge of 
the spectral content of the field being measured, allows a predetermined 
survey accuracy to be attained. Briefly, the two-dimensional (or in 
some applications, three-dimensional) spectral content estimates are 
used in algorithms which prescribe preferred track spacing, sampling in- 
tervals and track orientation. The method has been successfully applied 
to gravity, magnetic and physical oceanographic surveys. The availabil- 
ity of an adequate spectral content model for submarine topography would 
make this technology available for bathymetric survey design. 

In the detection of anomalous features in a field, the “normal” 
background variability must be removed by filtering to aid detection. 
This method has been applied successfully in magnetic anomaly detection 
and theoretically could be applied to bathymetric anomaly detection. 
McDonald, Katz and Faas (1966) applied this concept to submarine detec- 


tion, specifically in the search for the nuclear submarine Thresher. 


There are many applications of sea-floor topography modelling to 
the geologic interpretation of seafloor processes. Deterministic models 
illustrate that sea-floor relief is often decidedly anisotropic. This 
investigation will provide insight into the systematic nature of such 
spatial patterns and at several scales. Hayes and Conolly (1972) demon- 
strated the value of spectral techniques in studying the complex topog- 
raphy of the Australian-Antarctic Discordance. Since many geological 
processes result in linear features of characteristic wavelengths, an 
accurate, frequency-dependent, stochastic model could provide the capa- 
bility of decomposing such features and delineating their geographic 
distributions. It is anticipated that many important but totally 
unforeseen relationships will be discovered in exploiting this higher 
frequency model in the same manner that numerous fundamental relation- 
ships were discovered through the generation of global bathymetric 


charts. 
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3. Previous Work 


Although a great deal of quantitative geomorphology has been done 
on terrestrial landscape, as well as on lunar and planetary landscapes, 
relatively little quantitative study has been done on sea-floor morphol- 
ogy. Only since the late 1950's, when acoustic sonar systems became 
commercially available, has it been possible to attempt such statistical 
studies of bathymetry. The continuing refinement of sonar and naviga- 
tional systems has given recent investigators even greater opportunities 
for success. Equally important has been the development of large dig- 
ital computers and efficient statistical algorithms such as the fast 
Fourier transform, which allow sophisticated statistical analyses on 
large volumes of data, which were impractical until recently. 

The work of previous investigators is somewhat sparse and inconsis- 
tent, and is reviewed only briefly here. More information is given in 
later sections as it becomes pertinent to various aspects of the study. 
Some of the earliest studies were done by Agapova (1965), who generated 
mean, variance, skewness, and kurtosis statistics from measurements of 
slopes of a transect of the mid-Atlantic ridge. Heezen and Holcombe 
(1965) were able to describe physiographic provinces over a large por- 
tion of the North Atlantic Ocean. After rejecting spectral and filter- 
ing techniques, these authors developed a method of comparing the spac- 
ing of adjacent peaks and troughs. In effect, the method calculates the 
average distribution of slopes without regard to spatial frequency. 

Krause and Menard (1965) studied depth distributions from 15 pro- 


files in the east Pacific Ocean and found them to be normally distrib- 


uted. In studying normalized autocorrelation functions derived from the 
same profiles, the authors found no regularity with respect to lag. The 
profiles were used for delineating provinces on the basis of height ver- 
sus width ratio of the abyssal hills. Smith et al. (1965) continued the 
work of Krause and Menard (1965) and concluded that the distribution of 
depths is Gaussian when observed in distinct wavebands of 2 to 16 nauti- 
cal miles. Larson and Spiess (1970) later used a deep-towed instrument 
package to study the distribution of slopes in a small area of the east- 
ern North Pacific. Krause, Grim and Menard (1973) generated simple 
cumulative frequency plots of slopes in two areas of the East Pacific 
Rise. They found a very consistent power law form to describe these 
Ai Pronto and concluded that marine geomorphology could be described 
beers a few parameters. 

: Neidell (1966) generated spectral estimates of bathymetric, mag- 
netic and gravity profiles from the Atlantic and Indian Oceans. All 
spectra showed a "red-noise” character, that is, a decrease of power 
with increasing spatial frequency. The comparison of the various 
spectra was shown to be a valuable tool in investigation of complex 
geophysical problems. McDonald and Katz (1969) in their study of the 
directional dependence of roughness developed a polar autocorrelation 
function. Hayes and Conolly (1972) used spectral analysis as an inter- 
pretive tool in an area south of Australia. Distinct linear trends were 
interpreted by projecting randomly oriented tracks into north-south and 
east-west profiles and investigating the consistency of the resulting 
spectra. 

Clay and Leong (1974) rigorously described the relationship between 


surface roughness and the coherence of acoustic reflections. The dura- 


tions of returned pulses from hull-mounted sonars were used to estimate 
RMS roughness of microtopography (<0.2 km). Histograms of spectral 
estimates (periodograms) were generated by hand and shown to map con- 
sistently in an area southwest of Spain. 

Bell (1975b), also analyzing data from the eastern North Pacific 
Ocean, used power spectral techniques to generate composite estimates 
from several sources. He discovered a functional relationship for power 
versus spatial frequency of the form ax (power law form) to be consis- 
tent over a large range of enon scales. In a later paper (Bell, 
1978), this relationship was shown to hold for an enlarged data base, 
which is also confined to the same geographic area. The importance of 
anisotropy was recognized, and an initial look at the aspect ratios of 
submarine features was presented. 

Berkson (1975) generated spectra from a variety of bottom types and 
attempted to fit these with a power law form. Although a wide variety 
of coefficients were calculated, the power law form seemed to be consis- 
tent over many types of topography. Akal and Hovem (1978) generated 
two-dimensional spectra of sea-floor roughness from two sets of stereo- 
pair bottom photographs and a contoured bathymetric chart. They noted a 
remarkable consistency of form in all three spectra. Matthews (1980) 
developed a deterministic approach to describe bottom roughness. The 
North Atlantic and North Pacific Oceans were divided into 30 x 30 nau- 
tical mile squares and the maximum relief calculated. These cells were 
then grouped by range of relief (0-1100 m, 1100-1900 m, >1900 m) and the 
results mapped. Recently, Naudin and Prud'homme (1980) quantitatively 


described bottom morphologies from several areas based on multibeam 


sonar data collected by the SEABEAM system. Most recently, Berkson and 
Matthews (1983) have extended the work of Berkson (1975) and included 


estimates of sediment-basaltic interface roughness. 
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4. Statistical Considerations 


In order to develop a useful model of sea-floor roughness, one must 
select, from a seemingly infinite variety of available methods, an 
approach which is both suitable and tractable. In addition to consider- 
ing the nature of the sea floor itself, one must also consider such 
aspects as data resolution, computer storage capacity, statistical 
validity, and compatibility with various applications. Often, the 
choice of a particular method involves trade-offs between several of 
these considerations. In the following sections, many of these funda- 
mental considerations are addressed, and an initial approach to devel- 


oping this particular model is presented. 
Statistical Measures of Roughness 


If one considers surface roughness to be the variability of heights 

(or depths), the realm of statistics offers a multitude of measures to 
describe the roughness of a surface. Perhaps the most fundamental sta- 
tistic available to describe roughness would be one of the standard 
measures of data dispersion, such as the root mean square, standard 
deviation, or variance. These measures have the advantage of producing 
a simple parameter to describe the variability of depths in a given sam- 
ple. The major disadvantage of such simple measures is that they do not 

| provide information for roughness in terms of wavelength, and the sta- 
tistic decivedudecends upon the sample spacing and length of sampling as 


well as the actual roughness of the surface. 
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To illustrate the importance of having control over frequency 
dependence, consider two examples. If a generally flat surface which 
contains a high frequency roughness component of wavelength 4, were sam- 
pled at spacing A or any integer multiple thereof, each sample would 
fall at precisely the same depth and yield a variance of zero. Mathe- 


matically, all values of X would be identical, therefore the mean 


= n 

».4 =+ ie 2a would be equal to all X‘s and therefore 
roa (Xq - x)? 

Var (Xx) = eal = 0 


Although this example presents an extreme case, it is obvious that to 
assure an accurate measure of the variance at a given frequency, the 
sample spacing must be less than that corresponding wavelength, and the 
sample length long enough to sample all portions of the cycle. 

A more important shortcoming of these standard dispersion measures 
occurs at the longer wavelengths. Consider a generally smooth but 
broadly sloping surface. Examples from the deep sea might be a conti- 
mental rise or an abyssal fan. Since these dispersion measures record 
the average variation of individual samples from the sample mean, it is 
obvious that a relatively long sample would span a greater range of 
depths and produce a larger dispersion statistic than a smaller sample 
located in the same data. In this case, the measure of roughness would 
depend largely on the sample length. 

Many acoustic models of bottom interaction use the more general 
dispersion measure of the root mean square (RMS) roughness. In these 
acoustic models, this value represents the RMS variation of depth about 


some predicted value, normally the smoothed bathymetry. 


RMS roughness = [ 


where 2S represents a predicted value of depth at point i. By measuring 
the roughness relative to a predicted depth, rather than a simple mean 
as in the case of the standard deviation, the effect of long wavelength 
slopes on sampling is reduced. This improved measure does not, however, 
provide control over the distribution of roughness with frequency. 
Since the reflection or scattering of an incident acoustic signal on a 
surface is dependent upon spatial frequency, this often used parameter 
appears inadequate. 

Another possible measure of roughness is termed the “roughness 
coefficient” by Bloomfield (1976), and has the form 
2 
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Because this measure (also referred to as the von Neumann ratio and the 
Durbin-Watson statistic) is normalized by the total sum of squares of 
the residuals, the dependence on sample length is minimized. However, 
because this measure also depends on the squared differences of adjacent 
points, it measures only the variability of the signal at a wavelength 
corresponding to the sampling interval. This roughness coefficient then 
is essentially a ratio of the high frequency variability of a signal to 
its long-term trend, an interesting statistic, but not adequate for a 
complete stochastic model. 

Several statistical functions exist which describe the sample vari- 


ability as a function of discrete data spacings or lags (see Chatfield, 
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1980). Perhaps the simplest function of this type is the mean cross-—- 
product of terms at given lags k 


CK = =k ee Xq - Xg—ee = EL(X_)(X4-K)] 


This function (called simply the autocorrelation function in electrical 
engineering literature; see Bracewell, 1978) is both unnormalized and 
uncentered (the mean is not removed). Because it depends on absolute 
magnitude values (in this case, the total water depth), this simple 
measure can be improved for the purposes of roughness modelling by 


removing the sample mean, which yields the autocovariance function 
v(k) = E[(Xq - X) (X4-~ - X)] 


It is obvious that at a lag of k=0, (y(0)) is simply the sample vari- 
ance. By normalizing the autocovariance by the sample variance y(0), 


one derives the normalized autocorrelation function 
p(k) = ¥(K)/\ (9) 


Notice the equivalence between the normalized autocorrelation function 
at lag k=l, (p(1)) and the roughness coefficient described previously. 
In comparing the roughness of two different samples, it is undesirable 
to normalize by the sample variance, therefore, the autocovariance func- 
tion provides the best measure of variability with frequency (centered, 


but unnormalized). 
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Fully describing the variability of a process with its autocovari- 
ance function requires values at all n lags. By taking the Fourier 
transform of either the autocovariance or the unnormalized autocorrela- 
tion function, the process can be expressed more concisely as a function 
of frequency. This measure is the well known power spectral density 
function and can be estimated directly from the data with Fourler trans- 
forms. Besides providing a concise expression for the roughness as a 
function of frequency, the power spectrum also has many useful proper- 
ties which are described in Chapter 6 of Bracewell (1978). One theorem 
of particular interest is the derivative theorem which states that if a 
function f(x) has the transform F(s), then its derivative f'(x) has the 
transform i27sF(s). In this application, given the power spectrum of 
depths as a measure of roughness, the distribution of slopes (first 
derivative of depth) can be directly calculated. 

An additional measure of bottom roughness, which is particularly 
favored by those interested in acoustic modelling of bottom interaction, 
is the two-point conditional probability distribution function 
P(hy yhg|ry,F2)- This function defines the probability of measuring two 
heights (h, and hj) given two distance vectors (ry and FQ). Two consid- 
erations make this approach intractable. First, the description of the 
function requires a large number of parameters to be retained. Sec- 
ondly, complete two-dimensional data are required to adequately generate 
the function. This type of survey data is rarely available and only in 
areas which have been surveyed using multibeam sonar. 

The power spectral density function appears to be the best choice 
of statistical measure for bottom roughness, and it will underlie the 


stochastic model generated in this study. For convenience, the ampli- 
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tude spectrum (square root of power spectrum) will be used. When 
properly normalized, the amplitude spectrum allows the amplitude of 
component sinusoids to be expressed in simple length units, which can be 
interpreted more easily than units of length-squared. The method of 
calculation will follow Davis (1974) with proper windowing, prewhiten- 
ing, etc. As will be discussed in later sections, the simple and 
consistent functional form of spectra of topography allows relatively 
easy description and manipulation of the model. Also, recent work by 
Brown (1982, 1983) has shown the value of working in the frequency 


domain in modelling scattering from rough surfaces. 


Validity of Measurement Over Large Areas 


In generating the variance, autocovariance, or power spectrum from 
a discrete sample, only one realization of an infinite number of possi- 
ble realizations within the population is observed. The degree to which 
this realization is valid over the entire population depends upon the 
degree of homogeneity of the process. The condition of spatial homoge- 
neity is generally known as “stationarity” and the population being 
described referred to as a “stationary process”. The term “process” is 
used in the statistical sense of the variation of data with either time 
or space. In the case of sea-floor topography, the depths vary as a 
function of space. 

Stationarity is normally defined in two ways (see Chatfield, 1980; 
Popoulis, 1962). A spatial series is said to be “strictly” (or "“first- 
order”) stationary if the joint distribution of the process does not 


depend upon position. This implies that the mean and variance do not 
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depend on position. The less restrictive definition of “weakly” (or 
“second-order") stationary processes requires the mean to be constant 
and the autocovariance function to depend only on lag, not on position. 
This second definition has somewhat greater application, however, it is 
still too restrictive for general use with a spatial process as varied 
as submarine topography. 

Consider a broadly sloping surface, such as an abyssal fan. The 
mean depth in this province would by definition vary with position, and 
therefore would not satisfy the first requirement for stationarity. Yet 
if the process is homogeneous in higher spatial frequencies, one might 
prefer to treat this province as a homogeneous area for modelling. In 
practice the existing deterministic model could be used to describe the 
low-frequency trend. The sample data could be high-pass filtered to 
remove the non-stationary trend before generation of a spectrum. 

The presence of a low-frequency trend in samples of geophysical 
data is quite typical. In almost any length sample of a natural proc- 
ess, there are frequency components present with wavelengths greater 
than the sample length. In most natural systems, there is a finite 
limit on the rate of change of the process, causing the longer wave- 
length components to be also of greater amplitude. This typical spec- 
trum of natural processes was termed a “red-noise” spectrum by Shapiro 
and Ward (1960), an analogy to the red color of low frequency visible 
light. 

Although the red=-noise spectrum is the usual form in natural sys- 
tems, Figure 4-1 illustrates schematically that the presence of non- 
stationary components (in this case the statistical mean) can occur at 


any frequency and is dependent upon the horizontal extent of the sample. 
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C 


Figure 4-1 


Schematic illustration of the scale dependence of station- 
arity. All traces are derived from the same function (B), 
however the inferred stationarity of the mean would be dif- 
ferent depending upon the scale of observation. Observed 
at the scale of frame A, the mean is obviously non-station- 
ary. However, the mean appears stationary at the larger 
scale of B or the smaller scale of C. In most natural sys- 
tems, the non-stationary components occur in the partially 
resolved low spatial frequencies. This difficulty can be 


handled analytically by defining homogeneous provinces on 


the basis of stationarity and prewhitening of the spectra. 
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In order to generate a stochastic model of sea-floor topography, we must 
ensure that the process is stationary within some limit and over some 
defined area, with respect to the statistical parameters being calcu- 
lated. Ideally, this is accomplished by identifying homogeneous prov- 
inces based upon these criteria. Davis (1974) developed such a “prov- 
ince picker” for use in geophysical survey design, and his method is 
illustrated in Figure 4-2. 

By actively defining provinces which are weakly stationary in the 
frequency band of interest, one improves the validity of the statistical 
measures generated within each province. In addition, by delineating 
province boundaries, one can also alleviate the problems associated with 
the least-squares or averaging nature of Fourier transforms. In gener- 
ating a power spectrum from a data set, the operator must select the 
length and location of the sample series to be transformed. The result- 
ing spectrum will reflect the average frequency composition over the 
sample. If the sample spanned two distinct statistical processes, the 
result would provide the average composition of the two provinces, and 
would accurately represent neither. By confining one's samples within 
the boundaries of a homogeneous province, one insures a valid and repre- 
sentative statistic. These concepts will be discussed in detail in 


Chapter 5. 
Functional Representation of Spectra 
One property which makes the Fourier transform, both continuous and 


discrete, such a powerful analytical tool is its ability to express a 


spatial process in the spatial frequency domain both exactly and com- 


Band-Pass Filter (A) =(B) 


Low - Pass Filter (C) =(D) 


Contoured (D) = (E) 


Figure 4-2 Schematic illustration of the method developed by Davis 
(1974) for delineating homogeneous roughness provinces 
within a given spatial waveband. The original signal (A) 
is band-pass filtered at a pre-selected wave band of inter- 
est to produce B. The now centered data of B is full-wave 
rectified by taking the absolute value of all terms to pro- 
duce C. This output is low-pass filtered (smoothed) to 
yield the continuous function of D. Finally, this function 
is contoured at some selected interval (not necessarily 
linear) to delineate the boundaries of provinces X, Y, and 
Z in frame E. 
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pletely. This complete correspondence between domains allows detailed 
analyses to be made on the spectra with the assurance that there is an 
exact analog in the spatial domain. In order to maintain this corres- 
pondence, it is necessary to retain the complete transform, (both ampli- 
tude and phase components) in the spatial frequency domain. In the case 
of discrete data and transforms, it would be necessary to retain all of 
the degrees of freedom present in the original data set. 

In the creation of a probabilistic model, it makes little sense to 
retain as much information in the model as was present in the original 
data. Presumably, one would prefer to use the original data as a deter- 
ministic model. Also, the purpose of the model is to describe the gen- 
eral high frequency structure of the sea floor, rather than to analyze 
in detail the individual frequency components. Finally, the restric 
tions of computer storage space require a limited number of parameters 
in the model. 

All of the above considerations argue strongly for a severe paring 
of information in the frequency domain model. The phase spectrum, which 
requires fully half of the information in the spatial frequency domain, 
defines the origin in space of all component sinusoids of the amplitude 
spectrum and adds very little insight into the general structure of the 
sea-floor. An analysis of bathymetric phase spectra presented in Chap- 
ter 5 will show that the phase is in fact randomly distributed. For the 
purpose of this model, no phase spectra will be retained, as none of the 
previously stated applications require phase information. 

All measured data contain a component of random noise. Remotely 
sensed data are especially susceptible to measurement noise, the partic- 


ular noise problems in measuring oceanic depths having been mentioned in 
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the introduction. The presence of noise in the spatial data manifests 
itself in two ways in spatial frequency spectra. Inaccuracies of verti- 
cal measurements in the spatial domain result in the presence of a hori- 
zontal “white-noise level” in all amplitude or power spectra. This 
problem will be treated in detail in the following section. Uncertainty 
in the location of features in the spatial domain results in the scat- 
tering of amplitude estimates about the true frequency spectrum. These 
distinctions in the sources of error are somewhat artificial since the 
vertical and horizontal uncertainties are interdependent. 

Figure 4-3 reproduces a typical amplitude spectrum of depths. The 
red-noise character of the distribution as well as the scattering of 
amplitude values is apparent. The spectrum was derived from data col- 
lected by the U.S. Naval Oceanographic Office using SASS (Sonar Array 
Subsystem), and represents the highest resolution bathymetric informa- 
tion currently available from a surface ship (see Glenn, 1970). Control 
over relative horizontal location (navigational accuracy) of soundings 
is especially good due to the use of large, stable surveying platforms 
(in this case, USNS Dutton) and SINS (Ship's Inertial Navigation 
System). The degree of scattering of amplitude estimates would pre- 
sumably be greater in less sophisticated systems. 

Several methods come to mind to smooth this somewhat noisy spec- 
trum. A simple moving average taken over the amplitude estimates in the 
frequency spectrum would smooth the data. However, information would be 
lost from the high and low frequency extremes of the spectrum, while the 
density of data in the intermediate frequencies would not be reduced. 
The use of spectral windows or lag windows could be used both to smooth 


the spectrum and decrease the data density. The use of data windows in 
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Figure 4-3 
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Illustration of a typical amplitude spectrum of sea floor 
topography. A profile of the data, collected by SASS on 
the Gorda Rise, is presented above. The spectrum shows 
clearly the power law form as a linear fit on log-log plot. 
Possible extrapolations of the trend are shown as dashed 
lines beyond the fundamental and Nyquist frequencies. 
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the smoothing of spectra is discussed in many texts (see for example 
Bloomfield, 1976) and will not be discussed in detail here. The spec- 
trum presented in Figure 4-3 was in fact generated using the method of 
Davis (1974) which utilizes prewhitening as well as low-pass filtering 
of the spectral estimates. 

Another method of smoothing this rough spectrum is through the use 
of regression techniques. Calculating a continuous mathematical func- 
tion to describe the distribution of amplitudes would produce a smooth 
representation of the spectrum while, depending upon the complexity of 
the function used to fit the data, greatly reducing the number of param- 
eters retained in the model. This least-squares representation will be 
used in this study. 

By describing the spectrum with a simple, continuous mathematical 
function, one can more easily take advantage of the many symmetry prop- 
erties of the Fourier transform described by Bracewell (1978). For 
example, Rayleigh's Theorem (or Parseval's Formula for discrete series) 
states that the integral of the power spectrum equals the integral of 


the squared modulus of the function, or 
fP \£(x)|? ax = f° |F(s)|? as 
=O —o 


This is equivalent to stating that the total energy in one domain is 
exactly equal to the total energy in the other. If one were interested 
in the total energy in a particular band of frequency (in order to study 
the bottom interaction of sound of a particular wavelength, for exam- 
ple), the high and low frequencies of the band pass would become the 


limits of integration of the power spectrum. With the spectrum repre- 
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sented as a simple function, this definite integral could be evaluated 
analytically to derive the RMS variability for a discrete waveband. 

Having decided to use functional representations as the basis of a 
stochastic model of sea-floor roughness, the selection of a suitable 
functional form for the spectra becomes crucial. In order to minimize 
the size of the model, the simplest functional form which is justified 
by the data should be the best. An examination of Figure 4-3 as well as 
many other spectra presented later, would suggest the use of a simple 
straight line fit to the data. In light of the scatter in this already 
smoothed data, no higher order functional form is justified. 

The normal form for fitting a straight line to data is the estima- 
tion of the coefficients a and b in the equation 


a a 


y=bxta 


Notice in Figure 4-3, however that the data are plotted versus logarith- 


mic scales. The regression equation would therefore be written, 
log A =b log s + log a 


where A = amplitude and s = frequency. This equation can be rewritten 


in terms of A as, 
Aza. sd 


This inferred relationship between amplitude and frequency is often 


termed a “power law” or “power function” relationship. Appropriate 
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regression techniques must be selected in order to assure a proper 
regression fit to the power law form. Of prime concern is whether the 
residuals are minimized in log-log space or linear-linear space. The 
methods used in this study with accompanying theory and computer soft- 
ware are presented as Appendix A. 

The power law form seems to represent the best model for describing 
sea-floor topography in the spatial frequency domain. Its simplicity 
permits the frequency structure of a sample profile to be described 
using only two parameters, a considerable reduction of the original, 
deterministic data. Extensive work by Benoit Mandelbrot (1982) has 
produced a theoretical basis for this consistent relationship, formu- 
lated in terms of fractal dimension, a parameter functionally related to 
coefficient b above (see Berry and Lewis, 1980). Bell (1975b) discov- 
ered the same power law relationship in data from the Pacific Ocean that 
including lower spatial frequencies than those of interest in this study 
(see Fig. 4-4). Notice that Migure 4-4 plots power spectral density, 
rather than amplitude as in Figure 4-3 and other example spectra. 
Because the vertical axis in both plots is logarithmic, this exponentia- 
tion appears graphically as a linear transformation. Mathematically, 
the squaring of amplitude represents a simple doubling of the slope, 
1.e., multiplication of the b term by a factor of two. 

Bell (1975b) finds a fairly consistent relationship at many size 
scales, with the slope of the log transformed power spectrum of 6 = 
=2. Although this value probably represents a good mean estimate, the 
examination of many spectra in this study will show an accountable var- 
fation in this value. Berkson (1975) also discovered a significant 


variation of regression coefficients for spectra generated from differ- 
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Figure 4-4 
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Spectral estimates collected by T. H. Bell (1975b) from 
various sources. The linear distribution of the data 
illustrates the power law form of submarine topography 
spectra for a wide range of scales. Notice that this spec- 
trum is plotted versus power (amplitude squared) which 
results in a doubling of the ordinate and the slope of an 
amplitude spectrum. Refer to Bell (1975b) for data 
sources. 
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ent bottom types. Analyses presented in Chapter 5 will link this 
“universal” relationship 


a a 
Aza.s 


to the absence of stationary provincing techniques. 


Extension of Spectra into High Frequencies 


The component of random noise present in all empirical data, 
includes uncertainties from many possible sources in the total meas- 
urement system. There are uncertainties associated with the measuring 
device itself, for example, the errors in the timing of the arrival of a 
sonar pulse. The interference of external sources, such as radio waves, 
may also affect the measuring instruments. The nature of the environ- 
ment between the detection device and the process being sampled may 
introduce error, such as the variability of sea water sound velocity in 
bathymetric surveying. The truncation of significant digits in the 
storage of digital data introduces a finite level of noise, often called 
“round-off error.” All of these sources combine to form a total noise 
level for a measured data set. 

With the exception of round-off error, which is calculable, the 
source of these random errors can not be decomposed. However, using 
spectral techniques the level of the total noise can be estimated. It is 
a well-known property (see Bracewell, 1978, Chapter 16, for a complete 


discussion) that the spectrum of a randomly generated signal varies 
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about a constant value. This makes intuitive sense, since one would not 
expect any particular frequencies to dominate a random series. Again in 
analogy to the spectrum of visible light, this random component of a 
signal is often called “white noise”, reflecting the equal contribution 
of all frequencies. 

In the same way that Rayleigh’s or Parseval's theorems can be used 
to relate energy in one domain to energy in the other, the level of 
noise in a spatial signal can be estimated from the amplitude of the 
noise level of an amplitude spectrum. The noise level in the spatial 
domain is normally expressed as a simple dispersion measure of the vari- 
ability of the data, in this case, the oceanic depths. The following 
formula relates the white-noise level of the power spectrum to the root 


mean square. 


RMS = ¥ P/na 


where n = number of data points in series; and P = the mean power level 
of the white noise. In the case of sonar systems, knowledge of this RMS 
level defines the resolving capability of the system for a given signal 
level. Using these simple spectral techniques, the resolving power of 
the various sounding systems in use today can be calculated and 
compared. 

The red-noise structure of natural systems has been mentioned and 
illustrated previously (see Fig. 4-3). The interaction of natural sig- 
nals with instrument noise takes the form of a decreasing red-noise 
spectrum of the signal “intersecting” the horizontal white-noise level. 


In the spatial domain perspective, lower frequency features tend to have 
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higher amplitude and therefore the signal in these frequencies is “vis- 
ible” above the RMS noise. The frequency at which the spectrum of the 
signal intersects the noise level represents the highest frequency that 
is being resolved in the system. 

Figure 4-5 illustrates these concepts on a spectrum of SASS data. 
Notice that Figure 4-5 shows an obvious noise level while the spectrum 
shown in Mgure 4-3 does not. Both sets of data were collected using 
the same sounding system within days of each other, and it is expected 
that the instrument noise in each is approximately equal. The differ- 
ence is therefore that the data in Figure 4-3 represent a rougher area 
in which the signal energy in the highest frequencies is sufficient to 
maintain the resolution of the signal above the noise. The noise level 
is present in both sample spectra, however, it was never intersected in 
the sample from higher energy sea floor. The ability of a sonar to 
resolve horizontal features depends not only on the horizontal resolving 
power limitations of the instrument (normally limited by the size of the 
“footprint”), but also by the amplitude of features present in the sea 
floor. 

In light of the limitations of noise in all measurement systems, a 
fundamental question in the development of a stochastic model of sea- 
floor roughness from widely available surface sonar, is how far into the 
high frequencies the model can be extended. One could argue that due to 
the persistent exponential form of the spectra noted in this study, as 
well as the work of Bell (1975b), that this functional form can simply 
be mathematically extrapolated into higher frequencies. Further justi- 
fication might be provided by the generation of spectra from deep-towed 


instrument packages, bottom photographs, and direct observations, to 
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Figure 4-5 Illustration of an amplitude spectrum of sea floor topog- 
raphy which has encountered the "white-noise" level in the 
data. The spectrum of random noise has the form of a hori- 
zontal line. This noise level, which was derived from data 
collected by SASS, corresponds to a random component in the 
spatial data with root mean square dispersion of 1.9 
fathoms. 
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provide spot checks of the high frequency structure. This approach will 


be pursued in Chapter 7. 


Effect of Linear Features 


In generating a Fourier transform of topography within a stationary 
province, a one-dimensional statistic is generated from a two- 
dimensional surface. Provided the surface is isotropic, that is that 
there is no significant directional dependence, statistics derived from 
a one-dimensional sample would be equally valid for any orientation. 
Even a cursory examination of a deep-sea bathymetric chart reveals 
clearly that the sea floor, at least in the lower spatial frequencies 
presented in a contour chart, is quite anisotropic. There is extensive 
evidence that bathymetric lineations aieotertat to some extent in the 
higher spatial frequency roughness of interest to this study. To com- 
pletely model sea-floor roughness, it is essential to account for ‘any 
major directional dependence of the topographic features. 

Figure 4-6 illustrates the effect of sampling a simple two- 
dimensional periodic function in differing directions. Notice that the 
true wavelength (A) of the feature is sampled only when the sampling is 
exactly perpendicular to strike (8 = 0°). Any oblique angle produces an 
apparent wavelength (A') which is greater than A. At ® = 90° (a sample 
taken parallel to strike), the feature is not expressed at all (\' =©@). 


The apparent wavelength is related to the true wavelength by 


At = [costo] . A 
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Notice that the amplitude of the feature is not affected, provided a 
full wave form is sampled. 

The effect of directional sampling in the spatial frequency domain 
can be calculated using the previous relationship in combination with 
the similarity theorem of Fourier transforms (see Bracewell, 1978, 
p.- 101). The similarity theorem states that given the transform pair 


f(x) > F(s), then 
f(ax) > |a|~! F(s/a) 


Applying the geometry for directional sampling of linear features, i.e., 


a= |cos 6| 
£(|cos 8| . x) D |cos ellen - F (s/cos 9) 


Because |cos 6| must always be less than one, the effect of oblique 
sampling in the frequency domain is to shift the amplitude peak to lower 
frequencies, narrow its width, and increase its amplitude. This is a 
result of the fact that a signal of equal height but lower frequency has 
greater power than its higher frequency counterpart. It is important to 
note that these theorems assume an infinite continuous signal. In ana- 
lyzing finite length signals, it is necessary to normalize the spectrum 
by the sample length, that is, divide each amplitude estimate by the 
number of values in the time series. This normalization preserves the 
true amplitude of the individual waveforms, and allows comparisons 
between samples of different length. The effect of anisotropy on sam- 


pling a surface with continuous spectra is developed in a later section. 


34 


Several authors have recognized the importance of anisotropy of the 
sea floor to power spectral studfes. Hayes and Conolly (1972) recog- 
nized distinct peaks in their spectral analysis of the bathymetry of the 
Australian-Antarctic Discordance. These groups were normalized to two 
linear trends by projecting the data to simulate north-south and east- 
west samples. Distinct trends were successfully identified, however, in 
longer wavelengths than the high spatial frequency band of interest in 
this study. Bell (1978) also recognized the importance of anisotropic 
features in his study of abyssal hills in the north Pacific Ocean. One 
result reported in his study of the aspect ratios of features in this 
province is that the degree of anisotropy tends to decrease in the 
higher spatial frequencies. Whether or not this is a true relationship 
or the result of resolving limitations will be discussed further in 
Chapter 6. 

Figure 4-7 shows the sample locations for two spectra from the lin- 
ear Mendocino Fracture Zone in the northeastern Pacific Ocean. Figure 
4-8 presents the profile and amplitude spectrum for line A-A', which was 
sampled perpendicular to strike, reflecting the long wavelength fracture 
zone. Figure 4-9 is the corresponding plots for line B-B', which paral- 
lels strike and is located in the zone of disturbance. The exponential 
form of both spectra is evident, however Figure 4-10, which shows a com- 
posite of the trend of both ‘spectra, illustrates the differences in 
slope (exponent) of the two spectra. Segment A-A’ contains more power 
in the low spatial frequencies, while segment B=-B' contains more power 
in the high spatial frequencies. 

Both spectra are valid representations of the spatial frequency 


distribution in this physiographic province, but neither spectrum indi- 
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Figure 4-7 Bathymetric chart from the Mendocino Fracture Zone of the 
eastern Pacific Ocean showing the location of the data used 
in Figures 4-8 and 4-9. Notice that both samples lie 
within the same bottom environment, but that A-A' jis sam- 
pled across the long axis of the fracture while B-B' jis 


sampled along the axis. 
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Figure 4-8 Profile and amplitude spectrum of SASS data collected along 
line A-A' in Figure 4-7. 
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Figure 4-9 Profile and amplitude spectrum of SASS data collected along 
line B-B' in Figure 4-7. 
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Figure 4-10 Composite of Figures 4-8 and 4-9 illustrating the anisot- 
ropy of the sea floor in the Mendocino Fracture zone. 
Although both spectra retain their power law form, there is 
an obvious difference in slope and intercept. Profile B-B' 
appears to contain more high spatial frequency energy, 
while profile A-A' contains more energy in lower spatial 
frequencies. 
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vidually describes the roughness in all directions. The need for a 
directionally dependent function is obvious. The two-dimensional 
Fourier transform might appear to be appropriate, since it describes the 
two-dimensional roughness of a surface. However, its calculation 
requires a complete two-dimensional grid of data values which is gener- 
ally unavailable. The double Fourier transform, well described in Davis 
(1973), is calculated from two orthogonal one-dimensional spectra. This 
method, especially without the retention of the phase spectra, can not 
unambiguously identify the orientation of trend. Both the two-dimen- 
sional and double Fourier transforms require a large two-dimensional 
matrix to be retained in the model. 

McDonald and Katz (1969) describe a method for estimating autocor- 
relation functions as a function of azimuth 9. A similar approach will 
be attempted here for use with amplitude spectra. By studying the azi- 
muthally dependent distribution of the coefficients derived from the 
power law regression, a and b, it is anticipated that some functional 
form or forms will be revealed to allow modelling of the entire process 
as a function of both spatial frequency (s) and true azimuth (8). If 


a=f,(®) and b=f,() then 
F(8,8) = £,(6) . s*b(°) 


With this functional form, input into the model of simple compass 
direction for a given location would return the unique spatial- 
frequency-dependent function coefficients for the amplitude spectrum in 
that particular direction. Evaluation of these 9 dependent functions 
should provide a simple measure of the degree and direction of bottom 


anisotropy. 
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The derivation of the functional form of this 8 dependence would 
best be investigated and verified in areas where complete two- 
dimensional measures of depth are available. Areas surveyed with multi- 
beam sonar systems such as SASS and SEABEAM are ideal for this purpose. 
If a functional form is determined to be consistent, discrete samples 
from randomly oriented tracklines could be fitted with this functional 
model and F(®,s) estimated. Such studies will be presented in Chapter 6 
and Appendix E. 

One complication that could arise due to anisotropic roughness is 
in the use of the province picker for the delineation of stationary 
provinces. In a highly lineated area, the power level of particular 
marrow frequency bands would be directionally dependent. The province 
picker, however, measures total energy (integrated power) rather than 
discrete power, and as mentioned previously, the peak shift due to 
oblique sampling of a linear feature does not affect the total energy 
measured. Unless a major spectral peak is shifted beyond the low cut- 
off frequency of the band pass used in the province picker, the results 
should not be adversely affected by directional sampling problems. Even 
this problem seems unlikely to arise as thus far in all spectra sampled, 
none have shown unusual low spatial frequency peaks such as are seen by 
Hayes and Conolly (1972), and which are probably unique to a few tec- 


tonic provinces such as the Antarctic-Australian Discordance. 
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Summary 


In summarizing the preceding sections, a tentative strategy can be 


formulated for approaching the stochastic modelling of sea-floor 


roughness. 


(1) Delineations of homogeneous provinces, using a province 
picker, such as that developed by Davis (1974), would insure 
some degree of stationarity, and therefore validity, for the 
statistics describing the area. As will be discussed in the 
following chapter, this province-picking algorithm must be 
based on the same statistical measure that underlies the 
model, that is, the frequency spectrum. 

(2) Generation of amplitude spectra from available data within the 
delineated provinces would follow. Sample lengths for spectra 
generation would be confined to within the province boundaries 
defined in stage 1. One-dimensional spectra would describe 
the roughness in several orientations to provide input for 
later modelling of topographic anisotropy. 

(3) Regression analyses would be performed on these amplitude 
spectra to generate the coefficients of the functional form 
chosen to represent the spectra. Preliminary indications are 
that this form will be one or several power law relationships, 
each of which require only two coefficients to describe. 

(4) Anisotropy of the bottom would be modelled by studying the 
variation of the coefficients a and 6 of the power law model 


as a function of direction 8. In areas where complete infor- 
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(5) 


mation is available in two dimensions, such as a SASS or 
SEABEAM survey area, it might be possible to verify through 
regression analysis a simple ® dependent function to describe 
a and b. These functional forms, £,(9)=a and £,,(8)=b, could 
then define the best model for estimating model coefficients 


in areas where only randomly oriented track data is available. 


Extension of the functional representations beyond the spatial 


frequency at which the surface sonar encountered its white- 
noise level, would be attempted by studying the spectra of 
deep-towed sonar and bottom photographs. It is anticipated 
that a general functional form, perhaps a simple extrapolation 
of the power law form, will be discovered. In the many areas 
where high resolution bathymetry is not available, this 
extrapolated function should provide the best available esti- 


mate of spatial frequency structure at very short wavelengths. 


(6) Incorporation of the model into existing data sets would be 


the final development stage. This roughness model would be 
calculated for 5° grid cells and merged with the existing 
gridded data sets being developed at the Navai Oceanographic 
Office. All grid cells within a homogeneous province would be 
represented by coefficients generated from data located any- 
where within that province. By integrating this stochastic 
model with the existing deterministic models of oceanic 
depths, an essentially complete description of the sea floor 
will be contained in a single data base. This data base, when 
combined with similar gridded models of sea water sound veloc- 


ity, sediment column sound velocity, sediment thickness, and 
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other geophysical data, would provide a comprehensive environ- 
mental data base for further acoustical, oceanographic, geo- 


logical, and geophysical modelling and interpretation. 


The following chapters generally follow the approach presented 
above, beginning with a more rigorous look at the delineation of sta- 
tionary provinces. The method used for generating valid amplitude 
spectra are fully described in Davis (1974), and are reviewed only 
briefly. The regression techniques used in the study are described in 
detail in Appendix A. The problems of anisotropy and extension of the 
model into high spatial frequencies are discussed in detail in later 
chapters. Throughout, interpretation of the results in terms of geo- 


logical processes are presented. 
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5. Delineation of Statistically Homogeneous Provinces 


General Philosophy 


The importance of defining statistically homogeneous provinces was 
described in the previous chapter. The validity of frequency spectra, 
and indeed nearly all statistical measures, requires a stationary sample 
space. Unfortunately, truly stationary phenomena are usually only 
available to theoreticians and experimentalists. The statistics of most 
natural phenomena vary in either time, space, or both. It is therefore 
essential in attempting to describe statistically non-stationary phenom- 
ena, to delineate areas in which the statistic being generated varies 
minimally and only within defined limits. In order to accomplish this 
preliminary procedure of “province picking,” it is necessary to design a 
detection algorithm which takes account of the phenomenon being des- 
cribed and the statistic being used. 

To illustrate this principle of matching the province detector to 
the statistic being generated, we begin with an elementary statistic. 
As an example, assume that it is necessary to describe the areal distri- 
bution of height of the people of Africa. Assume for this discussion 
that the desired significance of the mean requires samples of at least 
10,000 individuals. One approach might be simply to divide the conti- 
nent into regular square areas and randomly select 10,000 heights from 
each area to generate a mean. The means so generated should indeed rep- 


resent the populations of these arbitrary squares. 
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To illustrate the effect of non-stationarity on the validity of the 
measured statistic, assume that in a particular sample square, the 
southern region is inhabited exclusively by Pygmies, averaging only 4 
feet tall, while the northern region is inhabited by Watusis, averaging 
7 feet tall. The me thod just described would predict a single popula- 
tion, averaging 5'6” height, inhabiting the area. In fact, very few of 
the individuals in the population are near this height and our statistic 
has failed to perform its intended function; to describe the areal dis- 
tribution of heights of the population. 

In order to confine our sampling to relatively homogeneous sample 
spaces, it is necessary to detect the boundary between independent pop- 
ulations prior to final sampling. The best method for accomplishing 
this: “province picking” is to measure the mean crudely with much smaller 
samples, say ten individuals or even one individual over smaller areas. 
Even these crude measures could be sufficient to define the large gra- 
dient in population mean across the boundary. Notice that the same sta- 
tistical measure (the mean) is used to describe the population and to 
define the province boundary. After the provinces are delineated, sam- 
pling within homogeneous provinces ensures statistical validity, at 
least to the degree that stationarity was confined in the province pick- 
ing procedure. An additional operational advantage of the province 
picking procedure is that very large areas of stationary means might be 
detected which would require only one random sample of 10,000 individ- 
uals to describe a large area, rather than conducting several repetitive 
samplings using the arbitrary grid technique. 

When one uses more advanced statistical measures to describe the 


earth, it becomes necessary to design more complex procedures for homo- 
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geneous province detection. The method used by Davis (1974) was 
described in Figure 4-2. In this application, the design of optimum 
survey spacing for marine gravity data collection, the statistic used 
was the total RMS energy in a particular spatial frequency band. The 
location of this band in frequency was dictated by later applications of 
the marine gravity data. Construction of a digital model of oceano- 
graphic sound speed requires delineating provinces in space and time 
based on statistics describing the shape of oceanographic profiles (T.M. 
Davis, personal communication, 1983). 

In order to delineate stationary provinces for the description of 
sea-floor roughness using frequency spectra, it becomes necessary to 
make a crude estimate of the amplitude spectrum discretely in the spa- 
tial domain. Recall that in transforming to the frequency domain, sta- 
tionarity has already been assumed, and therefore Fourler transform 
techniques are not appropriate. The method used in this study takes 
advantage of the relationship between band-limited energy in the spatial 
and frequency domains (Parseval's Formula), and the inferred power law 
form of amplitude spectra of topographic surfaces. Just as amplitude 
spectra represent the amplitude of component sinusoids at discrete fre- 
quencies, an equivalent estimate can be made in the spatial domain by 
band-pass filtering the frequency of interest and evaluating its ampli- 
tude. While the frequency domain estimate represents the least squares 
average amplitude over the entire leneen of sample, the amplitude: can be 
estimated discretely in the spatial domain using the Hilbert Transform. 
This is very similar to the method developed by Davis (1974) for a sin- 


gle frequency band. 
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In order to estimate the full spectrum, it is necessary to evaluate 
the amplitude at several frequencies, spanning the range of the desired 
spectrum. This is accomplished by convolving a bank of band pass fil- 
ters, centered at different frequencies, with the data and evaluating 
the amplitude of the band-limited signals discretely. Knowing the 
“power law” functional form of the spectrum in advance, one can fit the 
several amplitude versus frequency estimates at discrete points in 
space, using the iterative regression technique described in Appendix A. 
The regression coefficients a and b, mow available at every point along 
the profile, are often highly variable and must be smoothed. Also, 
because the two parameters are statistically correlated, it is prefer- 
able to use the exponent of frequency (b) and total band-limited RMS as 
detection parameters. Just as the presence of white noise at high fre- 
quencies must not be included in the regression analysis of the ampli- 
tude spectrum, amplitude estimates at the noise level in the spatial 
domain are also ignored. The method is described in detail in Appendix 
B, along with the results of various performance tests on signals of 


known properties used to calibrate the sensitivity of the detector. 


Generation of Amplitude Spectra 


Having delineated statistically homogeneous segments of data on the 
basis of their estimated frequency spectra, the next step is to generate 
amplitude spectra from these segments. Were the spatial domain esti- 
mates adequate, it would not be necessary to generate the spectra at 
all. However, due largely to instabilities in the slope (b) parameter, 


those estimates are not adequate and true FFI's must be run to estimate 
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the model parameters. Also, the Fourier transform method has the addi- 
tional advantage of estimating the amplitude at many more frequencies 
than the ten bands arbitrarily selected for the spatial domain 
algorithn. 

Since the proper generation of spectra is the basis for the entire 
model, great care has been taken to ensure that the best estimate of the 
true amplitude spectrum are obtained. The techniques used are described 
in detail by Davis (1974) and will only be reviewed here. The computer 
software used in this study was modified from programs provided by T.M. 
Davis and is presented as Appendix D. 

In using a finite length sample to represent an infinite series, 
the observer has in effect multiplied the infinite series by another 
infinite series consisting of zeros beyond the sample and ones at all 
sample locations. The multiplication of this so-called “boxcar” func- 
tion in the spatial domain, causes the true transform of the signal to 
be convolved with the boxcar's transform, a sinc function, in the fre- 
quency domain (see Bracewell, 1965). The presence in the frequency 
domain of side lobes on the sinc function, causes energy to be “leaked” 
into adjacent frequencies during convolution. Because of the red-noise 
character of spectra of sea-floor topography, this “leakage” tends to 
transfer energy artificially from lower to higher frequency. 

Although the use of tapered windows rather than boxcar sampling 
tends to reduce leakage, the preferred technique uses the method of pre- 
whitening. Tapered windows have a sinc function transform with eaduee’ 
sidelobes and a broadened mainlobe, which reduce spectral leakage at the 
expense of spectral resolution. In prewhitening, a specially designed 


high-pass filter is convolved with the signal, modifying it such that 
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its spectrum appears flat, or white, rather than red. When this pre- 
whitened signal is then passed through the Fourler transform, there is 
no preferential transfer of energy in either frequency direction. To 
obtain the “corrected” spectrum which approximates that of the true 
(Anfinite) signal, the prewhitened spectrum is divided by the impulse 
response of the prewhitening filter. This operation is equivalent to 
deconvolving the filter in the spatial domain. 

The importance of proper prewhitening can not be overstated. Leak- 
age of energy into high frequencies would cause a consistent underesti- 
mation of the magnitude of b (spectral slope), and degrade the ability 
of the model to estimate high frequency roughness (i.e., overestimation 
of amplitude at high frequencies). Figure 5-1 illustrates prewhitening 
by showing a raw spectrum, prewhitened spectrum, and corrected spectrum 


on one plot. Further examples can be found in Davis (1974). 


Physical Interpretation of Spectral Model Parameters 


Before examining the distribution of the spectral roughness model in 
selected study areas, it is worthwhile to discuss the physical meaning 
of the model parameters a and 6. The proportionality constant a in the 


expression 


a 


Azae° 3? 


where A = amplitude 


s = spatial frequency 
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Figure 5-1 Illustration of the importance of prewhitening of amplitude 


spectra. The raw spectrum is the result obtained by simply 
performing a spectral analysis on raw data. The prewhitened 
spectrum is the result of. performing the same analysis on 
data which has been convolved with a special high-pass fil- 
ter. The corrected spectrum results from the quotient of 
the prewhitened spectrum and the frequency spectrum of the 
high-pass filter, and represents an estimate of the spectrum 
of the corresponding infinite signal from which the sample 
data was derived. 
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represents a simple scaling factor for roughness. For a given frequency 
(s) and exponent (d), the amplitude (A) is proportioned to a. Due to 
the method of calculation of this expression, the actual value of a cor- 
responds to the amplitude of the component sinusoid with a wavelength of 
one kilometer and is usually expressed in meters or kilometers. This 
particular normalization was selected because the one kilometer wave- 
length falls within the sampling of most surface sonar data. For exam- 
ple, the required 0.5 km sample rate would be obtained with a 1 minute 
sonar ping rate on a ship traveling 30 km/hr (or 16 knots). The value 
of a does not necessarily correspond to any particular features in the 
signal, but only to the amplitude of the component sinusoid. 

The interpretation of the exponential parameter (b) is somewhat 
less intuitive. For the case of b= 0, the amplitude of all component 
frequencies is constant and equal to a. This ig the well-knowm “white 
noise” associated with random series such as instrument noise. Such a 
value for b would customarily be interpreted as instrument noise in any 
spectra from sea-floor topography. Values of b > O imply that ampli- 
tudes increase at shorter wavelengths, a condition that has never been 
observed in bathymetric data. What is consistently observed is the case 
where b < 0, the previously mentioned “red-noise” spectrum, in which 
amplitudes of component sinusoids increase with decreasing spatial fre- 
quency (longer wavelength). This indicates simply that broader features 
have greater height. 


An interesting special case occurs when b=-1. The expression 
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becomes 


or 


or 


A/A =a 


Simply stated, for b = -1, the ratio of amplitude to wavelength (or 
height to width) is constant and equal to a at all scales. This 
condition was termed “self-similarity” by Mandlebrot (1982) and corre- 
sponds to a fractal dimension of D= 1.5. 

One might visualize the special case b = -l as a signal which 
appears identical at all scales of observation. In another sense, the 
signal appears equally “rough” at all scales, the magnitude of roughness 
being prescribed by the magnitude of a. In cases where -l < b < 0, the 
ratio of height to width tends to decrease at longer wavelengths, 
although the absolute amplitude does increase. Such signals appear 
rougher at high frequencies. The converse case of b < - 1, implies that 
the ratio of height to width increases at longer wavelengths, and there- 
fore the signal appears smoother at high frequencies. Figure 5-2 sum- 


marizes these relationships. 
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Figure 5-2 Relationship of spectral slope parameter (b) to aspect ratio 
(A/A) of sinusoid at different frequencies. For b= -1, the 
aspect ratio remains constant in all frequencies. For b<-l, 
the aspect ratio increases at lower frequencies. For b>-1, 
the aspect ratio increases at higher frequencies. 
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Although the physical interpretation of the model parameters a and 
b is clear, an important question remains as to the geological signifi- 
cance of these terms. Why do certain areas of the sea floor have a par- 
ticular representative spectrum, and why do all spectra seem to show 
such a consistent power law form over large ranges of spatial frequency? 
An obvious hypothesis is that the spectral form reflects the unique 
interaction of the relief-forming processes and the materials being 
affected. For example, the formation of new sea-floor crust at oceanic 
ridge crests affects the relief of the new sea floor at all spatial fre- 
quencies. If the relief-forming process is uniform over some geographic 
region and interval of geologic time, there is no reason to suppose a 
change in the statistics of the surface being constructed, although its 
deterministic shape might change. Conversely, if there is a change in 
the reldief-forming process (such as the spreading rate) or material 
(perhaps a change in the properties of the magma source), it is likely 
that the resulting relief would also be affected. 

Many geological environments represent a composite of several 
relief-forming processes (tectonic, sedimentary, erosional) and several 
types of material. Such composite reliefs should result in an amplitude 
spectrum reflecting the composite spectra of these several processes and 
materials. If each style of relief is dominant over a different spatial 
frequency band, and each component spectrum conforms to the power law 
functional form observed in one-component cases, the composite spectrum 
should appear as a set of straight line segments on a plot of log ampli- 
tude versus log frequency. Examples of such composite spectra will be 


shown in a later section. 
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It is difficult to prove a direct relationship between statistical 
relief and combined process and material, but an interesting insight can 
be gained from a study of a highly variable environment, sedimentary 
microtopography. In 1981, Mark Wimbush of the University of Rhode 
Island deployed stereo camera on a structure on the upper continental 
rise northeast of Cape Hatteras. Stereo-pair bottom photographs were 
taken at an interval of twenty-seven days. Time-lapse photography 
showed that between the dates of these stereo-pair photographs, the fine 
scale sea-floor relief beneath the structure was altered both by biolog- 
ical activity and episodic bottom current events. 

Two microrelief maps were generated from the stereo-pair images and 
these are illustrated in Figures 5-3 and 5-4. Transects of heights were 
taken at 0.5 cm intervals (labelled A, B, C, D) across the surface. 
Amplitude spectra all showed the power law form found in spectra at 
lower frequencies, and in addition the spectral parameters showed no 
significant differences in spite of the gross change in the surface (see 
Figure 5-5). This implies that the two surfaces merely represent two 
realizations of the same statistical process. In terms of frequency 
domain analysis, only the phase spectrum is altered by the redistribu- 
tion of features, not the amplitude spectrum. This simple experiment 
does not unequivocally prove a causal relationship between statistical 


relief and process, but it does provide an encouraging result. 


The Phase Spectrum 


To reconstruct a profile or surface from its frequency domain repre- 


sentation, it is not sufficient to model only the amplitude of each com- 
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Figure 5-3 Contour representation of a bottom stereo-pair photograph 
collected on August 16, 1981. Dotted lines represent tran- 
sects used for spectral model generation. 
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Figure 5-4 Contour representation of a bottom stereo-pair photograph 
collected at the same location as that shown in Figure 5-3, 
but 27 days earlier, on July 20, 1981. The spectral models 
generated along indicated transects were not significantly 
different from those generated from transects of Figure 5-3. 
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Figure 5-5 Amplitude spectrum derived from profile A illustrated in 
Figure 5-3. Regression lines represent model spectra from 
profiles A,B,C, and D. Differences between these spectra 
are within estimation error, indicating no significant 
variation in time for the microtopography at this location 
nor significant anisotropy within each sample. 
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ponent frequency. One must, in addition, define the position of each 
component sinusoid relative to some geographic origin. The location of 
each sinusoid in space is expressed by its phase relative to this geo- 
graphic origin, and the composite of all component frequencies with 
their corresponding phases represents the phase spectrum. Since sines 
and cosines are trigonometric functions, phase is normally expressed as 
an angle between - 180° and 180°. 

Results from this study show that within statistically homogeneous 
provinces, the amplitude spectrum can be consistently modelled with a 
single or several power law functions. Although there is some varia- 
bility of the measured amplitude around the simplified model, the calcu- 
lated parameters remain consistent over often large geographic areas. 
However, any two sample profiles are not necessarily identical or even 
statistically correlated. The differences in the spatial domain man- 
ifestations of identical amplitude spectra can only be due to differ- 
ences in the phase spectra. 

Figure 5-6 illustrates a typical phase spectrum and the statistical 
distribution of its phase angles, derived from a single bathymetric pro- 
file. Several profiles were examined, which represented a variety of 
geographic locations and geological environments. The variability of 
phase angle with increasing frequency appeared to be random. A simple 
one-sample runs test was performed on several phase spectra, and all 
proved to be randomly ordered to within 95% confidence limits. The runs 
test is a non-parametric method (Siegle, 1956), meaning that no proba- 
bility distribution of the population is assumed. Examination of the 
distribution of phase angles indicates a uniform statistical distribu- 


tion; that is, a random distribution in which all phase angles are 
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Figure 5-6 Typical phase spectrum from a bathymetric profile collected 
in the Cascadia Basin. The upper diagram plots phase angles 
versus corresponding spatial frequencies. The histogram 
(below) shows the distribution of phase angles by 10-degree 
class intervals. Statistical testing indicates the series 
to be uniformly distributed random noise. 
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equally likely to occur for any frequency component. A Kolmogorov- 
Smirnov goodness-of-fit test was performed on several phase spectra, and 
mone were significantly (95%) different from the uniform distribution 
(Siegle, 1956). 

The spatial consistency of the amplitude spectrum and the uniformly 
distributed random nature of the phase spectrum of topography indicate 
that the differences in bathymetric surfaces within statistically homo- 
geneous provinces simply represent multiple realizations of the same 
statistical process. The observed changes in the microtopography 
recorded in Figures 5-3 and 5-4 can be modelled by combining two differ- 
ent random phase spectra of the same distribution with the known ampli- 
tude spectrum. With the functional representation of the amplitude 
spectrum for an area, a “typical” profile or surface can be produced by 
generating a uniformly distributed set of random numbers to represent 
the phase spectrum, and performing an inverse Fourier transformation to 
the space domain. 

The concept of representing the sea floor as a deterministic surface 
combined with stochastic variability was introduced in an earlier sec- 
tion. In that discussion, the deterministic components appeared as a 
smoothed, long wavelength surface which was combined with a higher fre- 
quency, stochastic roughness component. By describing the higher fre- 
quency components with a spectral representation, we can now visualize 
the amplitude spectrum as being determined (by modelling) and the phase 


spectrum as being a purely random (stochastic) process. 
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Creation of the Model and Interpretation 


Having developed the algorithms for defining quasi-stationary prov- 
inces and generating valid amplitude spectra, these methods were applied 
to an area off the Oregon coast. The area (42°N-45°N, 130°W-124°W) was 
selected due to data availability and the variety of geologic environ- 
ments represented within this relatively small area. The area includes 
the continental margin (shelf and slope), Astoria deep-sea fan, Tufts 
Abyssal Plain, Gorda Rise spreading center, Blanco Fracture Zone, the 
Cascadia Channel and numerous seamounts. All spectral estimates were 
generated from data collected on the SASS multibeam sonar system by the 
U.S. Naval Oceanographic Office. Only center beam depths were used in 
this portion of the analysis, in order to simplify processing. 

Figure 5-7 compiles the results of the combined province-picking 
and spectra-generating procedures. The areas delineated by the various 
shading patterns represent stationary provinces with similar ranges of 
the a statistic; that is, the amplitude of the component sinusoid at a 
wavelength of one kilometer. In cases of coincident values at crossing 
lines, a simple average was taken, ignoring for this analysis the 
effects of anisotropy. In many cases, provinces separated in the prov- 
ince-picking procedure became recombined in the final spectral model, 
indicating that the provincing algorithm used is more stringent than the 
levels selected for presentation. 

Within each of the larger provinces, the spectral slope parameter 
(b) estimates were averaged and these values shown within the provinces. 
The standard deviation of all estimates within any province was less 


than 0.1 in all but one case shown. There is one province in which the 
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Figure 5-7 Distribution of roughness in the vicinity of the Gorda Rise, Northwest Pacific Ocean. 
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slope parameter is different in two areas within a single shading pat- 
tern (bottom center of chart). Many roughness province boundaries coin- 
cide with obvious physiographic province boundaries. In these cases, 
the bathymetry was used to trace the province boundaries between sample 
spectra. In many other cases, no boundary was obvious in the bathymetry 
and such tracing was not possible; that is, what appears as an abrupt 
boundary may represent simply an arbitrary contour of a continuous gra- 
dient. Notice also that the illustrated bathymetry from Chase et al. 
(1981) was based on a totally different data set than that used for 
spectral model generation, and therefore some provinces appear in the 
model which are apparently unsupported independently in the bathymetry. 

The gross distribution of the roughness statistic a corresponds 
fairly well with what one would expect intuitively. The roughest areas 
(a > 2.5 m) are located at the Gorda Rise crest, Blanco Fracture Zone 
and over most seamounts, in particular the President Jackson Seamounts. 
The least rough areas correspond to the sedimentary provinces of the 
Tufts Abyssal Plain, Cascadia Basin, the Astoria deep-sea fan, and the 
continental shelf. The Cascadia Channel appears as an intermediate 
roughness province which can be traced very easily through the Blanco 
Fracture Zone and onto the Tufts Abyssal Plain. The Astoria Channel, 
located to the east of the Cascadia Channel, is too narrow (<8 km) for 
this analysis and thus does not appear as a separate province. 

Within this gross distribution of roughness, some more subtle pat- 
terns can be identified. The continental margin between the shelf and 
abyssal plain appears to be banded with the topography becoming gener- 
ally rougher down slope. This particular continental slope represents a 


slowly converging margin between the North American Plate and the Juan 
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de Fuca-Gorda Plate. Carson (1977) and Barnard (1978) have both exam- 
ined this convergent margin in areas to the north off the coast of 
Washington, and both present seismic cross-sections of this area. Pre- 
sumably, similar processes are at work along the Oregon margin. Barnard 
(1978) infers a change of compression rate from 2.3 em/year before 0.5 
mybp to a present rate of 0.7 cm/year. Deformation of Cascadia Basin 
sediments is progressing westward, making the deepest areas of the mar- 
gin also the most recently deformed. This westward progression of 
deformation process is expressed in the bathymetry as a downslope 
increase in surface roughness. Barnard (1978) classifies the slope ter- 
rain into an upper slope extending to a depth of about 1500 m, and an 
accretionary “borderlands” complex of en-enchelon, anticlinal ridges. 
Between these ridges are sediment-filled basins. These physiographic 
divisions, derived by qualitative observation of the geological struc- 
ture of the region, correspond closely to the roughness model generated 
by quantitative methods. 

The Gorda Rise is one of the more active areas of the world sea 
floor, and a corresponding complexity is evident in the derived pattern 
of bottom roughness provinces. Atwater and Mudie (1973) reviewed the 
tectonic history of the area. Additional history of spreading rate and 
spreading direction changes can be found in Elvers et al. (1973). Much 
of this tectonic history may be peripheral to this study because the 
sea floor affected is now buried beneath the Tufts Abyssal Plain and 
Gorda Deep-Sea Fan sediments. 

One interesting feature of the bottom roughness chart presented in 
Figure 5-7 is the very rough ridge crest which terminates abruptly on 


either flank. The full width of the feature is about 25 km. Were 


ridge- forming processes constant through time and ridge crest relief 
“frozen” into the topography, the same roughness would be expected to 
persist, at least in long wavelengths, on older sea floor. The other 
major process affecting the ridge flanks, sedimentation, would be 
expected to affect the short wavelengths initially (due to their lower 
amplitude), and form a smooth transition zone, not the abrupt boundary 
observed in the roughness pattern. Magnetic data from the area indicate 
that this zone falls within the Bruhnes-Matayama boundary and must 
therefore represent crust younger than 0.7 my. Recalling that Barnard 
(1978) found evidence of a slowing of compression rate on the continen- 
tal margin at a time younger than 0.5 mybp in areas to the north, it is 
possible that this roughness feature reflects the same change in proc- 
ess, most likely a slowing of spreading rate. Future examination of 
other ridge axes should reveal whether this pattern is unique to the 
Gorda Rise or present under other tectonic conditions. 

Perhaps equally interesting is the abrupt termination of this ridge 
crest at latitude 42°20'°N. The roughness values drop (as supported by 
three tracklines) two and three roughness levels at the feature's termi- 
nus. This disruption of the ridge crest falls along a trend which 
encompasses President Jackson Seamounts and other seamounts to the 
northwest, and a major (900 m) bathymetric deep to the southeast. The 
break in the ridge crest trend also appears in the bathymetric chart. 

Hey (1977) developed a “propagating rift” model to describe the 
plate geometrics and magnetic anomaly pattern found on the Juan de Fuca 
Ridge by the Pioneer survey. According to this model, the growing 
spreading center propagates along strike as the dying spreading center 


becomes inactive and is added to one of the rigid plates. As this proc- 
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ess continues through geologic time, a V-pattern of fossil spreading 
centers forms a “propagator wake”, as is illustrated in Figure 5-8. 
This same pattern is found in the isosynchronous magnetic anomaly pat- 
tern. Since this model was proposed, similar processes have been 
observed on other spreading centers, in particular the Cocos-Nazca 
spreading center (Searle and Hey, 1983). 

A propagating ridge crest model offers one explanation for the 
abrupt termination of the Gorda Rise crest shown in Figure 5-7. In 
order to test this hypothesis, a magnetic anomaly map of the Gorda Rise 
was constructed and is illustrated in Figure 5-9. The chart is based on 
original magnetic anomaly data collected by the U.S. Naval Oceanographic 
Office. The V-pattern associated with the “propagator wake” is evident, 
extending to the east and nae tmaselin the ridge crest at 42°N. The 
direction of the V-pattern indicates that the ridge crest to the north 
is propagating toward the south at the expense of the southern portion 
of the Gorda Rise crest. The geometry of the schematic model by Hey 
(1977) for the. Juan de Fuca Ridge (Mgure 5-8), represents a nearly per- 
fect analog to the magnetic anomaly pattern of Gorda Rise (Figure 5-9). 
It would appear that the abrupt termination of the high roughness zone 
shown in Mgure 5-7, is due to its association with the tip of a prop- 
agating rift system. 

Another feature of interest in the distribution of the a parameters 
in Figure 5-7 is the existence of east-west trends of selected roughness 
provinces on the ridge flank. One might have expected roughness prov- 
inces to align themselves sub-parallel to the ridge strike, perhaps 
reflecting changes in processes through time being felt along the length 


of the ridge axis. This is the case with the very rough central valley 
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Figure 5-8 Schematic illustration (from Hey, 1977) of the pattern of 
isochronous seafloor resulting from a southward propagating 
rift. Double lines represent active spreading centers, 
dashed lines are fossil spreading centers, heavy line is 
active transform fault, and dotted lines are associated 
fracture zones. Diagonal trend lines indicate the so-called 
“propagator wake”. 
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Figure 5-9 Magnetic anomaly chart of the Gorda Rise area. Positive 
anomalies are shown in black; negative anomalies are shown 
in white. The northeast-southwest trending Gorda Rise is 
obvious as is the northwest-southeast trend of the Blanco 
Fracture Zone near the northern limit of the chart. Compare 
the inferred “propagator wake" (indicated by dashed lines) 
to the theoretical model of a propagating rift illustrated 
in Figure 5-8. 
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portion. The east-west trend of roughness indicates instead relief- 
forming processes which act relatively constantly through time, but 
quite variably along the ridge strike. Francheteau and Ballard (1982) 
describe the along-strike variability of processes along the East 
Pacific Rise, and relate the changes to distance from the ridge crest/ 
fracture zone intersection. The change in process is expressed by the 
relative importance of fluid lava flows and pillow lava flows. These 
petrologic changes are in turn associated with the elevation of the rift 
valley along the ridge crest; topographic highs are associated with high 
ratios of fluid lavas and topographic lows associated with pillow lavas. 
Indeed the shallowest portion of the Gorda Ridge segment is located near 
latitude 42°45'N, which shows relatively low roughness values as one 
would expect from the sheet-like flow of fluid lavas. Ridge flank areas 
become rougher adjacent to deeper axial valley segments, which might 
reflect the rougher surface of pillow lavas. Although other explana- 
tions for the trend of roughness on the ridge flanks could be put forth, 
the distributions found in this study are consistent, although not nec- 
essarily typical. Extension of the model to more thoroughly investi- 
gated ridge crests should shed light on this particular hypothesis. 

The patterns apparent in the distribution of a are not as evident 
when one examines the distribution of the spectral slope (b). It is 
clear that the “universal” value of -1 for the spectral slope inferred 
by Bell (1975b) and others is not supported by this modelling effort. 
The reason for this discrepancy is not readily apparent, however the 
attention given to defining stationary sample space in this study does 


represent one major difference in method. In order to test this hypoth- 
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esis, amplitude spectra were generated for several profiles in the Gorda 
Rise area which spanned multiple stationary provinces. The effect of 
these extensions was to reduce the degree of statistical homogeneity of 
each profile. The computed spectral slope (b) parameters for the 
resulting amplitude spectra converge consistently to the value b = -1.0 
as longer profiles are tested. Figure 5-10 illustrates one such long 
profile extending nearly 500 km from the Oregon coast. The profiles 
analyzed by Bell (1975b) were often much longer. No formal statistical 
argument for this convergence to b = -1.0 will be attempted here, how- 
ever, it would appear that the concatenation of multiple profiles of 
differing spectral characteristics results in a profile which resembles 
a random walk. The need to define statistically homogeneous sample 
spaces before generating statistics is clearly demonstrated. 

Most slope values shown in Figure 5-7 cluster about -1.5 with the 
exception of the smooth ridge axis segment south of 42°20'N. The two 
large sedimentary provinces are represented by two values for spectral 
slope. Figure 5-11 illustrates a typical amplitude spectrum from these 
sedimentary provinces. The spectrum clearly separates into two distinct 
straight-line segments of different slope. The average values of these 
distinct slopes for all profiles is given in the corresponding boxes 
(see Figure 5-7). If the hypothesis is accepted that the characteristic 
spectral slope of amplitude spectra of sea-floor topography represents a 
dominant relief-forming process, these spectra should represent areas 
where two processes are at work, affecting the relief in different spa- 
tial frequency bands. It is likely that the higher frequency band rang- 
ing from A = 2.5 km and with 6 = -.6, represents the sedimentary regime 


in these areas. The lower frequency process with b < - 1.4, 1s less 
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Figure 5-10 Amplitude spectrum of a long bathymetric profile, which 
encompasses several statistically homogeneous provinces. 
The inclusion of non-stationary segments into the analyses 
results in a spectral slope parameter which approaches b= 
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Typical amplitude spectrum of profiles collected on the 
Tufts Abyssal Plain. Spectrum shows two distinct 
straight-line segments which intersect at a spatial 
frequency of ~0.4 cycles/km. The longer wavelength 
portion is described by = -1.847, a = 0.0294 m. The 
shorter wavelength segment is described by 6 = -0.541, 4 = 
0.113 m. Profiles from the Cascadia Basin province have 
similar spectra, however the model parameters are slightly 
different. 
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obvious, but may well represent an underlying tectonic effect due per- 
haps' to the compressional nature of the area overprinting a long wave- 
length relief on the smooth sediment. 

This general correspondence of b < -1 in tectonic provinces, and -1l 
< b < 0 for sedimentary provinces is found in many areas outside this 
study area. The reason for this general relationship can only be specu- 
lated upon. Recall from the previous section that b < -1 requires that 
the ratio of height to width of component features, increase at longer 
wavelengths. If one envisions tectonic relief forming processes, the 
entire morphology is constructed and then erosional and sedimentary 
processes begin smoothing small features first, progressing to con- 
stantly larger scales. In sedimentary processes, one can envision a 
smooth layer of sediment being affected by bottom current interaction 
for example. This constructive process begins at the highest spatial 
frequencies and progresses to larger scales. This is in agreement with 
the relationship of -1 < b< 0, in which the ratio of height to width of 
component features increases at shorter wavelengths. 

The study area illustrated in Figure 5-7 was selected because of 
its rich geological diversity. As such, the distribution of roughness 
provinces is correspondingly complex. To allow a comparison with a more 
tectonically stable area of the world ocean floor, a large number of 
tracklines off the United States east coast were analyzed to compare a 
passive continental margin. Although the trackline density was not ade- 
quate for a complete chart to be drawn, one interesting relationship was 
discovered. An extremely large area of the continental margin, compris- 
ing most of the continental rise, was found to have in common a very 


distinct amplitude spectrum. Figure 5-12 illustrates a typical ampli- 
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Figure 5-12 Typical amplitude spectrum of profiles collected on the 
Continental Rise, east coast U.S. All fourteen profiles 
examined in this province (shown in Figure 5-13) show 
nearly identical patterns. The two linear segments, in 
this example intersect at ~0.3 cycles/km., with 6 = 
-1.813, 4 = 0.0142 m. for the longer wavelength model, and 

= -.524, a = 0.0645 m. for the shorter wavelength 
model. 


tude. spectrum from this large province, and shows the same two-process 
mature of the spectrum in Figure 5-ll. Im fact, this spectrum is almost 
identical to the measured spectra from the Tufts Abyssal Plain and the 
Cascadia Basin. 

Figure 5-13 shows the location of profiles with this characteristic 
spectrum. Some of these profiles are over 200 nom long, and the profiles 
are distributed over an area more than 1000 nm in extent. The total 
standard deviation for the spectral parameters in all fourteen profiles 
was only 0.15 for the slope (b) parameters and 0.02 meters for the 
intercept (a) parameters. These values apply to both the lower fre- 
quency and higher frequency line segments. The area may well extend 
even further northeast or southwest. Profiles are only broken by sea- 
mounts or deep=sea channels associated with submarine canyons. As such, 
one can see that the areal extent of a roughness province with a given 
level of allowed non-stationarity, can wary, econ hundreds of thousands 


of square miles to the slopes of a single seamount. 
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6. Anisotropy of Surfaces 


The necessity of directionally treating anisotropic surfaces was 
introduced in Chapter 4. Any statistic generated from a one-dimensional 
profile of a two-dimensional surface is only valid in all directions if 
that surface is isotropic. This chapter examines the importance of 
anisotropy in detail. A simplified theoretical model of the effect of 
anisotropy on the frequency spectra of directionally sampled profiles is 
formulated and tested. Such spectra are generated for two areas of the 
sea floor where complete, two-dimensional bathymetric data are available 
from multibeam sonar. An identical study is performed on data from 
side-sean sonar. These results are then compared to the theoretical 
model. Next, the possibility of estimating such two-dimensional func- 


tions from randomly oriented bathymetric profile data is discussed. 


Theoretical Model 


Before examining the effect of anisotropy on measured frequency 
spectra derived from actual bathymetric profiles, it is instructive to 
examine a very simple theoretical model of this effect. Several of the 
concepts introduced in Chapter 4 will be utilized. There it was shown 
that the effect of one-dimensionally sampling a sinusoid which has been 
extended to two dimensions (see Figure 4-6), is to stretch the true 


wavelength as, 


AY = |cos™!6| Oi 
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where A" = apparent wavelength 
4 = true wavelength 
@ = angle of sampling (0° = perpendicular 


to linear trend) 


This relationship can then be combined with the similarity theorem 


of Fourier Transforms to yield the transform pair 


£ (|cos 6] * x) D> |cose|~! + F (s/cose) 


Recall that these relationships were formulated for a single component 
wave form extended to two dimensions. Actual sea-floor topography has a 


spectrum which {fs continuous and conforms to a power law functional form 


Amzae gb 


To extend the model to the continuous case, one can envision gener- 
ating a topographic profile (or other signal with continuous power law 
form, such as a random walk model), and extending all points on the pro- 
file to the second dimension. This surface is then sampled in various 


directions and the spectrum of each profile evaluated for a and b above. 
Combining the above relationships yields 
£(|cos 0| * x) D |eos e|7! + a = (s/cos 0)” 


or equivalently 


£ (|cos 0|~! = x) D |cos 8| + a+ (s * cos @)> 
or 


F (s, 0) = |cos 0| + a * (s * cos 6)? 


where F(s, 8) is the Fourier transform of £(|cose|~1 ° x) expressed 
in terms of 6 and s. The coefficient a can now be expressed as a func- 


tion of 6 as 
a(6) = |cos0| * a 
which can be peneralized to 
a(@) = |cos(6-8,)| * a 


where 6, {s the azimuth perpendicular to the linear trend. Notice 
that for © = 6, that is a profile generated perpendicular to trend, 
|cos(0-6,)| = 1 and F(s,0) = a * 8, the original function. For 6=0, = 
£90°, |cos( - 94) | = 0° and F(s,9) = 0, that is, for profiles sampled 
parallel to strike, the series is a constant (normalized to zero) and 
therefore contains no energy. The parameter b is independent of 0. 

Figure 6-1 illustrates the hypothetical “plane wave” surface that is 
used in this simple model. As a test of the above theory, radial tran- 
sects of the surface were generated and the resulting profiles input to 
the standard Fourier transform routines used throughout this study. 


Figure 6-2 depicts the sampling patterns used to generate the profiles. 
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Graphic representation of an elementary theoretical model 
for anisotropic surfaces. The surface is created by gener- 
es a ue oi walk (with theoretical amplitude spectrum of 

= ) in the x (90° azimuth) direction. These values 
are met extended to the second dimension, creating a 


lineated surface with strike = 0°. Viewpoint is from the 
southwest. 
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37 TRANSECTS 
EVERY 5° 


Figure 6-2 Radial sampling pattern used to generate one-dimensional 
amplitude spectra from various azimuths on a surface. 
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Figure 6-3 plots the value of a (“intercept”) and b ("slope of spec- 
trum”) as a function of azimuth. The resulting spectral estimates are 
fitted with the above model. Approaching azimuths of 0° and 180°, the 
profiles nearly parallel the linear trend and the results are unreli- 
able, due to the small number of depths available for analysis resulting 
in a very narrow frequency band width used in the model regression. 

The coefficients b(6) do not show any relationship to the trend, as 
predicted by theory. This does not imply that directional dependence in 
b(6) is never found in spectra from sea-floor profiles. If a surface 
had two or more distinct signals (perhaps due to differing relief- 
forming processes) superimposed, cyclical behavior in b(@) would be pos- 
sible. For example, envision a hypothetical surface composed of an iso- 
tropic two-dimensional signal with spectra slope b,(@), overlain by a 
simple linear trend (like the one described above) with spectral slope 
bg(®). Perpendicular to trend, b(8) would be some combination of b,(6) 
and bg(®) depending on their relative amplitudes. Parallel to trend, 
b(@) would equal just b, (8), in this example, since the linear trend 
with slope bg(9) is constant in the direction parallel to strike. The 
example spectra from the Mendocino Fracture Zone presented in Chapter 4 
illustrate this effect. Appendix E examines a series of artificially 
generated surfaces and their spectral characteristics. 

In the results shown in Figure 6-3, the parameter a(@) shows the 
expected relationship to cos(6-6) 0, = 90°, as predicted by theory. 
This model cosine function can be used to parameterize the surface 
anisotropy. The model sinusoid is generated by an iterative regression 


technique (see Appendix C.2) which determines the following equation 
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Distribution of spectral parameters versus azimuth of samp- 
ling for theoretical surface shown in Figure 6-1. The 
upper series represents the slope of the spectrum in log- 
log space and varies randomly around the theoretical value 
of -1 as predicted by theory. Near 0° and 180°, the esti- 
mates become unstable due to poor sampling. The lower ser- 
jes represents the intercept of the spectrum in log-log 
space (or the coefficient of frequency) and agrees well 
with the sinusoid model predicted by theory. Notice the 
maximum intercept, and therefore total spectral energy, 
corresponds to a sample taken perpendicular to strike 
(azimuth = 90°). 
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a(9) =u +v° cos(2 ° (9-8,)) 


the three term regression technique yields estimates for u, v, and Qo. 
These terms have definite physical meaning. u represents the simple 
mean roughness level of the surface, that is the mean a(@) of the signal 
sampled in all directions. This could be visualized as the “isotropic” 
component of the surface. v determines the amplitude of the sinusoidal 
component of the regression model and represents a measure of the degree 
of anisotropy of the surface. 80 estimates the normal to the true azi- 
muth of the linear trend. Frequency is not estimated since the perio- 
dicity of 1 cycle/180° is known. 

Unfortunately, it is not possible to decompose more than one linear 
trend in a surface using this method. Envision a surface consisting of 
two linear trends of differing orientation (8 a and 6g), “anisotropy” 
levels (va and vg), and “isotropy” levels (ug and ug)- The surface, 


being a simple linear combination of the two component trends can be 


expressed as 
a(®@) = (ug + ug) + va © cos(2 * (6-0,)) + vp ° cos(2 * 8&-6g)) 


In this example, u, and ug are both presumably equal to zero for “per- 
fect” linear trends. However, even in non-perfect cases in which some 
energies are available parallel to strike, the u components are linearly 
combined and can not be differentiated. The “anisotropy” components 
also combine linearly to yield another sinusoid whose amplitude and 


phase are dependent on the relative amplitudes and phases of the orig- 


inal sinusoidal components. Appendix C presents a geometric proof of 
this relationship. 

In the case where two or more linear trends are present in a sur- 
face, the form of a(@) will show a simple sinusoid with phase and ampli- 
tude which can not be uniquely decomposed into component sinusoids. If 
an estimate of the azimuth and amplitude of one of the trends could be 
produced independently, this component could be removed and the remain- 
ing component sinusoids analyzed. It should be emphasized that the 
inability to decompose component trends in no way invalidates the model, 
it simply complicates the interpretation of the model statistics in 
terms of formation processes. 

Multiple linear trends can only be decomposed in cases where the 
trends are sufficiently band-limited to appear as distinct peaks in the 
frequency spectrum. This approach allows decomposed trends to be 
uniquely identified by spectra generated in two orthogonal directions, 
as was shown by Hayes and Conolly (1972). Their work clearly shows such 
trends in the large scale topography of the Antarctic-Australian 
Discordance. In examining a great many spectra of small-scale (i <8 km) 
topography during the present study, no significant spectral peaks have 
been observed, even in topography which is highly lineated at the larger 
scales. This result may be due in part to the presentation of these 


spectra in log-log form. 
Comparison of Theoretical Functional Forms with Multibeam Sonar Data 


The previous section developed a simple theoretical model of the 


effect of linear trends on frequency spectra from profiles sampled at 
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varying azimuths on an anisotropic surface. Figures 4-7-10 illustrated 
the effect of linear features on two bathymetric profiles collected at 
near right angles. To test the validity of the theoretical model to 
actual bathymetry, ie is necessary to sample regularly around the con- 
pass at a single location on the sea floor. This is only possible for 
areas which have complete areal bathymetric coverage of high spatial 
resolution. 

The most practical instruments available for obtaining such data 
are the multibeam sonar systems. These systems, which provide a “swath” 
of discrete soundings on a line perpendicular to the ship's track, allow 
complete coverage of an area. By conducting surveys in which the paral- 
lel survey tracks are spaced so that the outer beams of adjacent tracks 
overlap or are nearly juxtaposed, a complete two-dimensional survey can 
be performed. Very few of such data sets are currently available. How- 
ever, two data sets from contrasting geologic environments were made 
available for this study (see Chapter 5, Section B). Both surveys were 
conducted by the U.S. Naval Oceanographic Office using the SASS multi- 
beam sonar system (Glenn, 1970). The recent acquisition of the academic 
SEABEAM systems should provide more full-coverage surveys in the future. 
New techniques for processing side-scan sonar data from the SEAMARC-1 
system allow similar two-dimensional analyses to be performed at smaller 
scales. 

Figure 6-4 presents the contoured bathymetry from the Gorda Rise 
area of the northeast Pacific Ocean. Contours were created automati- 
cally and only appear where supported by multibeam soundings. The 
coverage is generally complete with the exception of a small area on the 


western margin of the chart and small gaps between swaths. This chart 
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Figure 6-4 Index maps showing the location of two-dimensional SASS 
bathymetry data used in study of azimuthal dependence of 
topographic spectra. The data are in the vicinity of the 
crest of the Gorda Rise in highly lineated topography. 
Central coordinates are 42°53.5'N, 126°40'W. Contour 
interval is 10 fathoms. 
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was created using a grid spacing of 0.25 minutes of longitude and lati- 
tude (~460m x ~300m), well above the resolving capability of the SASS 
system. 

The area shown in Figure 6-4 was selected to test the effect of 
anisotropy on one-dimensional amplitude spectra. The data set is 
located in the vicinity of the ridge crest, which was identified as a 
quasi-stationary province by the methods described in Chapter 5. The 
lineation of the topography is obvious on the index chart and trends 
approximately N25°E. 

Figure 6-5 represents graphically the test area. Although the grid 
spacing used for the illustration is 0.1 minutes of latitude and longi- 
tude, the profiles were generated from a grid with spacing of 0.05 min- 
utes (~100m), which approaches the resolving limit of the SASS system 
for these water.depths and noise level conditions. Again the sampling 
pattern shown in Figure 6-2 was used to generate radial profiles of 256 
points. No ensemble averaging was used. 

The spectral parameters a(®) and b(6) are plotted versus azimuth in 
Figure 6-6. The results agree closely with the theoretical model devel- 
oped in the previous section. Notice first that the parameter b(6), 
plotted above in these diagrams and labelled “Slope of Spectrum", shows 
no systematic fluctuation with azimuth, as predicted by theory. Note 
also that the mean slope is -1.24, well below the -1.0 slope of the ran- 
dom walk (Markov process) model. The large variability of this paran- 
eter could be reduced by ensemble averaging of several spectral esti- 
mates, created by offsetting the center of the sampling pattern. The RMS 
variability would be reduced by 1/ Y N , where N is the number of esti- 


mates (see Chapter 7). 


SASS Bathymetry — Gorda Rise 


Figure 6-5 Graphic representation of SASS bathymetry data projected 
onto an evenly-spaced grid. Illustration uses 128 x 128 
points spaced at 0.1 minutes of latitude and longitude 
without cartographic projection. Fourier analysis was per- 
formed on 256 points from a 0.05 minute grid. Measured 
strike of the lineations is approximately 25° azimuth. 
Viewpoint is from the northeast. 
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Figure 6-6 Distribution of spectral parameters versus azimuth for 


Gorda Rise spreading center bathymetry shown in Figure 6-4 
and 6-5. Spectral slope parameter (above) shows no appar- 
ent functional relationship to azimith as predicted by 
theory (notice that the mean slope is -1.24, well below 
-1.0, which would correspond to a Markov process). The 
intercept parameter (below) clearly shows the effect of 
seafloor anisotropy and generally conforms to the sinu- 
soidal model. Model parameters are as follows: mean ampli- 
tude = 1.68 m, amplitude of sinusoid = 0.51 m, azimuth of 
maximum energy = 115° (perpendicular to observed strike). 
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The results for a(8) also conform reasonably well to the model pre- 
diction. The values, plotted below and labelled “Intercept”, represent 
the amplitude (in meters) of the Fourier component at a wavelength of 1 
kilometer. The mean intercept is 1.68 meters, indicating a relatively 
rough topography, and corresponding to the “isotropic” term u in the 
model. The anisotropy of the surface is obvious from the large ampli- 
tude (0.51 meters) of the model sinusoid (the v term). The maximum 
value of the model occurs at an azimuth of a, = 115°. This corresponds 
to the normal to the linear trend (which was measured as 25° in the 
full-coverage chart) as expected. These values fully parameterize the 
effect of surface anisotropy on the frequency domain description. 

Although the functional models for b(®) and a(®) appear to be of 
the proper form, there are some obvious variations in the measured par- 
ameters from the model. Im order to give an intuitive impression for 
the degree that the generalized model departs from the actual measured 
spectrum, Figure 6-7 illustrates a “worst case” example in which the 
modelled a(@) differs from the observed a(@) by ~0.5 m. The selected 
profile is from 9 = 115° (see Figure 6-6). Plotted with the measured 
spectrum is the regression line derived from the functional model, 
rather than the least-square fit usually shown. The model-derived spec- 
trum provides an excellent representation of the spectrum and appears to 
fall well within the estimation noise of the spectrum, even in this 
“worst case” example. 

The second study area is shown in Figure 6-8. This data set repre- 
sents a totally different style of sea-floor topography due to the geo- 
logic environment, which is controlled by sedimentological processes 


rather than the tectonic setting of the Gorda Rise. The broad trend of 
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Figure 6-7 Example amplitude spectrum from the Gorda Rise crest sampled 
at azimuth 115 degrees. The straight line segment through 
the spectrum represents the model prediction, rather than a 
least squares fitted line. As seen in Figure 6-6, this 
particular azimuth represents a large deviation of the 
fitted parameters from the model sinusoid. The model line 
still appears to fall well within estimation noise, even for 
this “worst case” example. 
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slope toward the south-southeast is very long wavelength and not 
included in the analysis. The small channels traversing the slope are 
due presumably to downslope transport of sediment and represent an 
apparent linear trend with wavelength of approximately 10-15 kilometers. 
This is at the low frequency limit of the present analysis, but might 
indicate such a trend in shorter wavelengths. Higure 6-9 graphically 
illustrates the data set (in this case gridded at 0.1 minutes of lat- 
{tude and longitude), and shows this linear trend due to down-slope 
processes. Coincidentally, the survey track was run quasi-parallel to 
this trend which could further complicate interpretation. 

The distribution of spectral parameters with azimuth are plotted in 
Figure 6-10, in the same format as the previous plots. Notice again 
that the parameter b(8) plotted above shows no functional relationship 
to azimuth 8, as predicted by theory. In this case, the mean value of 
the b(®8)'s is -1.05. The intercept parameter a(@) reflects the rel- 
atively smooth, isotropic nature of this sample of the sea floor. In 
this case, the mean amplitude of 0.097 meters is less than 62 of the 
same parameter in the Gorda Rise area. The “anisotropy” term, vin this 
ease is estimated at 0.0063 meters, only 1% of the value for the Gorda 
Rise. With these almost isotropic conditions, the estimate of the azi- 
muth of maximum energy cannot be made with any fidelity. 

The final test area represents an intermediate level of both gen- 
eral roughness and degree of anisotropy. The data were collected by the 
SEAMARC-1 side-scan sonar system, which is a deep-towed instrument 
developed by W.B.F. Ryan of Lamont=Doherty Geological Observatory. The 
vehicle is towed at approximately 500 m above the sea floor and collects 


data from side-scan sonar and from down-looking sonar. The depth of the 
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Figure 6-9 Graphic representation of SASS bathymetry data projected 
onto a 128 x 128 point grid. Grid spacing is 0.1 minutes 
of latitude and longitude and is presented without cartog- 
raphic projection. Visible in the surface are spikes assoc- 
jated with sonar processing noise and smooth areas where 
surface was interpolated between tracks. Viewpoint is from 
the southwest. 
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Figure 6-10 Distribution of spectral parameters versus azimuth for con- 


tinental rise bathymetry shown in Figures 6-8 and 6-9. 
Spectral slope (above) shows no apparent functional rela- 
tionship to azimuth and has a mean slope of -1.05. The 
intercept parameter (below) reflects the relatively smooth, 
isotropic nature of the surface. Model parameters are as 
follows: mean amplitude = 0.097 m, amplitude of sinuscid = 
0.0063 m, azimuth of maximum energy = 120°. Due to the low 
ive) of anisotropy, the azimuth direction is not signif- 
cant. 


vehicle is continuously measured. Farre and Ryan (1984) have developed 
a method of combining these sources of information into a detailed con- 
tour chart of bathymetry. These contours, when evaluated for depth on 
an evenly spaced grid, comprise the data set used in this study. 

Figure 6-11 illustrates a contour chart of the study area, which 
encompasses the Carteret Canyon, continental slope, and upper continen- 
tal rise. The smaller area outlined was used to produce the spectral 
parameters illustrated in Figure 6-12. The data grid uses a spacing of 
only 5 meters and represents much higher resolution than the surface- 
derived sonar data in the other two examples. Unfortunately, due to 
round-off errors producing a white-noise level of 0.3 meters in the data 
(only whole meter depths were retained), most spectra were only above 
noise to the 100 m wavelength band, similar to the resolution of SASS. 

Figure 6-12 illustrates the results of the azimuthal dependence of 
spectra study. The slope parameter (b) shows a mean of -1.88, lower 
than the other two study areas. The data also indicate what might be an 
example of azimuthal dependence of b, such as that observed on the 
Mendocino Fracture Zone and discussed in Chapter 4 and Appendix E. The 
intercept (a) parameters are u = .87 Mm, v = .23 m and relative azimuth 
(8,) = 85°, which represent a level of roughness and anisotropy inter- 
mediate to the previous two examples. Because the original survey was 
collected at an azimuth of 140°, the true azimuth of anisotropy is 45°. 

In conclusion, the simple theoretical model of the effect of linear 
trends on frequency domain statistics, can adequately describe anisot- 
ropy in three test areas. The parameter b(6) appears to be independent 
of azimuth in at least two areas, although quite noisy. The sinusoidal 


construction of the intercept parameter as 
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Figure 6-12 


SEAMARC-1—CONTINENTAL SLOPE 
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Distribution of spectral parameters versus azimuth for con- 
tinental slope/submarine canyon bathymetry shown in Figure 
6-11. Spectral slope (above) has a mean value of -1.88, 
although it may also contain a cyclical component not noted 
in other examples. The intercept parameters (below) repre- 
sent an intermediate level of both mean roughness and 
anisotropy between those of the Gorda Rise and Continental 
Rise examples. Model parameters are as follows: mean 
amplitude = 0.87 m, amplitude of sinusoid = 0.23 m, azimuth 
of maximum energy = 85°. Due to the direction of sampling, 
true azimuth on the earth corresponds to 45°. 
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a(@) =u +v¥° cos(2 * (0-89)) 


allows the background, or isotropic, roughness, as well as the degree of 
anisotropy and its trend to be quantified, and therefore compared in 
different areas. 

If one assumes this simple model of anisotropy to be true, at least 
in some cases, an interesting insight into the effect of scale on 
anisotropy can be seen in the mathematics. The assumption of a constant 
value of b(@) at all azimuths can be envisioned as a family of lines in 
log-log space of constant slope whose levels vary with azimuth. The 
orthogonal azimuths which represent the extremes of anisotropy, would 
maintain a constant spacing in amplitude at all frequencies in log-log 
Space. That is, if at a given frequency the amplitude in one direction 
were twice that of the normal azimuth, that relationship would remain 
constant at all frequencies. 

In examining bathymetry and other geological data, it often appears 
that anisotropy decreases at shorter wavelengths. Bell (1975) reached 
that conclusion in studying the aspect ratio of shapes of seamounts of 
different sizes. Much of the validity of this statement depends upon 
how anisotropy is defined. As stated above, if a relative doubling of 
amplitude in the orthogonal direction occurs at one scale, this same 
doubling should occur at all scales. However, if one considers the 
absolute difference in amplitudes in orthogonal directions, at long 
wavelengths the difference between perhaps one and two meters of ampli- 
tude represents one meter of difference, while at very short wavelengths 
this same relationship might appear as the difference between one and 


two centimeters. Although the proportional relationship of amplitude 
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(doubling) is constant, the absolute difference decreases exponentially, 
depending upon the value of b. Since the ability to resolve amplitude 
is limited, the ability to resolve anisotropy at small scales is also 
limited and could lead to an erroneous conclusion concerning anisotropy 
at high frequencies. These relationships do not apply in cases where b 


varies regularly with azimuth, such as those discussed in Appendix E. 


Estimation of Two-Dimensional Spectra from Randomly 


Oriented Ship Track 


An alternative method for describing the topographic roughness of a 
surface over all azimuths is with the two-dimensional amplitude spec- 
trum. The method involved is quite similar to that used in the one- 
dimensional case, but prewhitening requires a circularly symmetric high- 
pass filter, and a two-dimensional Fast Fourier Transform algorithm is 
used. Perhaps most important to the practical use of this method is the 
requirement for a complete two-dimensional array of data (depths) as 
input. 

As mentioned previously, complete areal bathymetric surveys are 
available in very few areas of the world ocean. To be practical, such 
surveys must use a multibeam sonar array such as SASS or SEABEAM, and 
tracks must be spaced so that adjacent swaths are juxtaposed. Of the 
areas presented in the previous section, the Gorda Rise survey shown in 
Figure 6-4 indicated the highest degree of anisotropy and was therefore 
selected for two-dimensional FFT analysis. 

Figure 6-13 illustrates the results of generating a two-dimensional 


amplitude spectrum from the gridded bathymetric data shown in Figure 
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Figure 6-13 Two-dimensional amplitude spectrum from the Gorda Rise area 
illustrated in Figure 6-5. Log-transformed amplitude esti- 
mates, computed via a two-dimensional Fourier transform, 
are represented by light contours drawn every 0.5 order of 
magnitude. The heavy lines represent amplitude estimates 
predicted by the four-parameter, azimuthally dependent 
model derived for the area. 
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6-5. Amplitude estimates appear as irregular contours. The amplitude 
values were log transformed before plotting, and these values are plot- 
ted at integer increments. Contours are drawn each .5 order of magni- 
tude of amplitude. The data set, and therefore its amplitude spectrum, 
is oriented with the columns parallel to longitude and rows parallel to 
latitude. 

Plotted with the spectrum in heavy lines is the two-dimensional 
spectrum as modelled from one-dimensional profiles by the method 
described in the previous section. Because the gridded data base used 
in the analysis was spaced evenly in latitude and longitude, the spec- 
trum as a function of spatial frequency is necessarily distorted. High 
frequency noise associated with the east-west oriented track lines 
appears as a smearing of the contours in the vertical and horizontal. 

The simple model spectrum expisines most of the variance in ampli- 
tude. The lineation of the topography with a strike of 6 = 25° can be 
easily seen as an elongation of the contours in the cross-strike (@ = 
115°) direction. The degree of anisotropy (v) term in the model deter- 
mines the elongation of the contours. The model appears to overestimate 
the amplitudes in the high frequencies slightly, which would indicate a 
slightly lower slope (b) value than that derived by the described 
method. An improved fit results if the value b= -1.5, the value 
derived for the ridge crest in Chapter 5, is used. Overall, however, 
the match of the true two-dimensional spectrum with the model is quite 
good. It is questionable whether the additional detail present in the 
true spectrum represents a valuable signal or simply additional noise. 

There are several advantages to the model proposed in this study 


(derived with respect to azimuth) over the two-dimensional FFT method 
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(constructed in cartesian coordinates). The proposed model requires 
only four parameters (b, u, v and ee) to describe the surface rough- 
ness. The two-dimensional spectrum in this case requires a 128 x 128 
array, or 16,384 parameters. In addition, the four parameters used in 
the model have physical meaning attached to them which may prove useful 
in comparisons of different areas. Also, the computer algorithms used 
to generate the model require far fewer calculations than the direct 
transform method. 

Perhaps the most relevant advantage in the azimuthal model con- 
struction is that the two dimensional nature of the surface can be esti- 
mated from randomly oriented ship tracks. Each profile yields a one- 
dimensional estimate of the amplitude spectrum at the azimuth of the 
ship's heading. Such estimates can be thought of as cross sections 
through the surface contoured in Figure 6-13. For example, the profile 
and amplitude spectrum shown in Figure 6-7 represent a cross section of 
the two-dimensional surface collected at azimuth NIL15°E. Given a suf- 
ficient number of such randomly oriented tracks over an adequate range 
of headings, the model can be constructed as described in this chapter. 
Until a great many more multibeam surveys have been collected, this 
method of estimating the anisotropy of bottom roughness will remain the 


only available method over most of the world oceans. 
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7. Prediction of High Frequency Roughness 


Having developed a spectral model of sea-floor topography based on 
measurements from surface ship sonar systems, the question remains 
whether this model can be extrapolated into spatial scales smaller than 
those resolved by the sounding system. This question is particularly 
important to underwater acoustic applications, where the acoustic fre- 
quencies of interest in a scattering problem do not necessarily corre- 
spond to the spatial frequencies sampled to generate the model. The 
concepts of measurement noise levels were introduced in Chapter 4. In 
this chapter, the effect of estimation errors on prediction will be 
examined, sources of measured high frequency bathymetry introduced, and 


a simple prediction test presented. 


Source of Error in Spectral Estimates 


The model parameters used to describe the amplitude spectrum of sea- 
floor topography are derived from regression estimates of spectra from 
profiles of noisy data collected in a generally non-stationary environ- 
ment. As such, there is estimation error in the model parameters from 
several sources, which will necessarily result in prediction errors as 
the model is extrapolated to high frequencies. Although many techniques 
are used to reduce these errors, some level of error will always remain. 
Unfortunately, due to the variability of data quality, variable track- 
line spacing, and the presence of some level of non-stationarity in a 


data province, it is not possible to quantify completely the estimation 
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error associated with the model. Only by comparing estimated values 
with values measured in high frequencies over many data sets and geo- 
logic environments, can a reasonable statistical base be assembled to 
assess the prediction capabilities of the model quantitatively. 

The effect of instrument noise (illustrated in Mgure 4-5) when 
encountered by the signal spectrum, has the effect of reducing the 
spatial frequency range of amplitude estimates available to the regres- 
sion analysis. Certain spectra examined in the course of this study 
varied from complete white-noise spectra to spectra showing only two or 
three amplitude estimates above the noise. Such spectra are of no use 
in model generation. The length of the spectrum available for regres- 
sion analysis depends as well on the length of data available to the 
FFT. Long profiles from large quasi-stationary provinces produce cor- 
respondingly long amplitude spectra and therefore more reliable model 
estimates. As mentioned in Appendix B, a minimum of one hundred points 
is allowed in a profile for analysis. The use of a higher minimum pro- 
file length, while improving regression estimates, would allow more non- 
stationarity in the profile and reduce the ability to resolve small 
provinces. 

No attempt has been made to quantify the level of non-stationarity 
(as defined for this study) by examining the mean variability of the 
spatial domain estimates used in province picking (see Appendix B). Due 
to the nature of the sea floor, certain provinces appear over thousands 
of square miles, while others fall smaller than the hundred point mini- 
mum required for analysis and must be combined. Such variability prob- 
ably precludes estimating the degree of non-stationarity with any accu- 


racy; we can only attempt to constrain the effect with the province 
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picking procedure. In general, however, one can treat estimates from 
large provinces of persistent geological processes as more reliable than 
those generated in a relatively small area. 

One source of error in the model can be quantified, and that is the 
residual error from the regression. Using standard statistical tech- 
niques (see, for example, Draper and Smith, 1981), the difference in 
individual amplitude estimates from the regression model estimates can 
be expressed as a root mean square. Errors in this study averaged €, = 
t .03 and E1og ae = .015m for a single spectrum. These errors reflect 
many of the errors associated with the model, although they cannot be 
decomposed into component sources. 

Another important factor determining the level of error in the 
model is the number of estimates used in generating a composite spec- 
trum. In the case of an anisotropic area, a variety of azimuths dis- 
tributed about the compass allows a better estimate to be made. In gen- 
eral, the ensembling of N time series composed of signal with noise, 
results in a decrease of noise (as RMS) of 1/ Vv N. In the case of our 
spectral model, the signals are the derived amplitude spectra along an 
azimuth. A simple method of improving the prediction capability and 
accuracy of this model is to ensemble-average the amplitude versus fre- 
quency estimates from several proximal and near-parallel tracks. The 
multibeam sonar provides exactly this capability and future developments 
should take advantage of it. 

Figure 7-1 illustrates this reduction of estimation error through 
ensemble-averaging of multibeam sonar derived spectra. For this exam- 
ple, sixteen parallel beams (profiles) from the SASS multibeam system 


were analyzed for the statistically homogeneous province of the Gorda 
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Figure 7-1 Results of ensembling the spectral estimates (a and b) from 
sixteen parallel profiles collected by a multibeam sonar 
system. Plotted points indicate the estimated parameters 
for each of the sixteen profiles and various ensembles of 
adjacent profile parameters. The calculated standard error 
(solid line) decreases as VN (dashed line), as predicted by 
standard statistical theory. 
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Rise crest. Each derived E-W profile was 17 km in length and spaced 
approximately 100 m apart. Amplitude spectra were generated for each 
profile, and rather than averaging each amplitude estimate, the derived 
regression parameters a and b were assembled for adjacent beams in 
groups of 2, 4, 8, and 16. Figure 7-1 plots for both a and b, the 
estimated parameters for each of the sixteen profiles and the various 
results of ensembling. The standard error of each set is plotted as a 
solid line, with the theoretical :!/VN relationship shown as a dashed 
line. The derived standard error for the final average of sixteen 
points is necessarily zero. Similar techniques could be applied in the 


province picking algorithm to improve reliability. 
Propagation of Error to High Frequency Estimates 


Although the errors associated with the model can not be determined 
with any accuracy, it is still instructive to examine how the errors 
affect the ability to predict amplitude in frequencies beyond the range 
of analysis. Due to the power law form of the model, 

A= a 0 ab where A = amplitude 
s = spatial frequency 


a a 


a, b = regression parameters 
the errors of estimate (error bars) are not linear. It is therefore 


somewhat easier to visualize the function in log A-log s space where the 


error bars are linear. 
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As previously stated, the parameter a represents the intercept (in 
meters) of the function with frequency log s = 0, or wavelength ’4 = 1 
km. The error associated with a (e,) appears as a constant vertical 
shift of the regression line in log-log space. It is in fact a multi- 
plicative factor in linear space and has the effect of multiplying or 
dividing the A value by ([antilog e,! at any frequency. Since €, is 
{ndependent of frequency, it is stable over large extrapolations. Since 
all spectra of topography are red noise, the absolute level of estima- 
tion error effectively decreases at higher frequencies (lower ampli- 
tudes). 

The spectral slope parameter b te not independent of frequency. In 
log-log space, the dimensionless bd appears as the slope of the linear 
regression line. Error associated with b (ey) at s = Q, causes an 
increasing prediction error at higher or lower frequencies. The rela- 
tionship in linear space is also multiplicative and depends on frequency 
as multiplying or dividing the A value by [antilog (|log 3| * €,)]. The 
total error of estimate requires linearly combining the two sources €. 
and €,, which translates into multiplying or dividing the value A(s) by 


[antilog (c€, + |log s| ° €,)]. An example of these calculations is 


included in the following section. 


Comparison of Surface Ship Sonar Results to Deep-Towed 


Sonar Results and Results from Bottom Photography 


To quantify accurately the ability of the spectral model derive 
from surface ship sonar systems to predict amplitudes at high spatial 


frequencies requires a large data base of small-scale bathymetry pro- 
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files. At present, such a large data base does not exist, at least for 
the meter-millimeter scales of bottom photography. One area of the 
world ocean (near 40°27'N and 62°20'W) has been extensively surveyed at 
small scales and this is the location of the High Energy Benthic 
Boundary Layer Experiment (HEBBLE). 

The HEBBLE area falls into the large roughness province of the East 
Coast Continental Rise which was described in Chapter 5 and illustrated 
in Figures 5-12 and 5-13. A spectral model was generated for this large 
area based on averaging spectral estimates from all profiles illustrated 
on Figure 5-13. As mentioned previously, the amplitude spectrum for 
this area consists of two segments, the lower frequency model (with a = 
0.0142 a, b = -1.813) extending to wavelengths of approximately A = 3 
km, and the higher frequency model (with a = 0.0602 m, b= -0.603) 
extending from A = 3 km to A = 200 a. 

Data collected by the Deep-Tow sonar system in the HEBBLE area were 
provided by Scripps Institute of Oceanography. The Deep-Tow, which col- 
lects profiles from a height of only 25-50 meters above the sea floor, 
is able to sample bathymetry at a horizontal sample spacing of 5 meters. 
The vehicle is positioned via a transponder navigation system which pro- 
vides relative location to the beacons every five minutes (or approxi- 
mately 270 meters) of track. All depths recorded by the system in this 
area appear to be above the instrument noise level. 

Digital height data from one stereo-pair bottom photograph col- 
lected in the HEBBLE area were provided by the University of Washington. 
The data set consists of six horizontal and six vertical transects sam- 
pled at 1 mm spacing. All derived spectra were virtually identical in 


their model characteristics. 
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Figure 7-2 presents a composite amplitude spectrum showing data 
from all three sources in the HEBBLE area. The SASS derived spectrum 
(from the nearest available trackline) spans wavelengths of 100 km to 
200 m. The model regression lines, derived from spectra generated 
throughout this large province, are shown as straight line segments. 
The Deep-Tow derived spectrum, spanning wavelengths of 1 km to 10 on, is 
very well predicted by the model regression line. The lower frequency 
estimates begin to diverge at A= 300 m., which may be due to positioning 
distortions from the transponder navigation system, which is interro- 
gated at 5 minute (or ~300 meter) intervals. Finally, the bottom photo- 
graph derived spectrum, spanning wavelengths of 25.6 cm to 2 mm, falls 
approximately .5 orders of magnitude below the regression line. 

The prediction residual of .5 orders of magnitude over 5 decades of 
frequency indicates some combination of errors in parameter estimates a 
and d. If all errors were in Ay it would be in error by .5 orders of 
magnitude. Similarly, the error in b would be 2.1, using the relation- 
ships described in the last section. The evaluation of several more 
bottom photograph-derived profiles from the same region would allow a 
more quantitative treatment of the model prediction error. It is pos- 
sible that other high frequency spectra may scatter around the predic- 
tion line, due perhaps to actual variability of the microrelief in the 
area rather than error in the prediction model. Visual inspection of 
the several photograph pairs taken in HEBBLE indicate that other areas 
are in fact rougher than the relatively featureless photo analyzed here 
(Arthur Nowell, personal communication, 1983). 

In any case, the prediction of millimeter scale topographic ampli- 


tude from a surface hull-mounted sonar system to within half an order of 
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Composite amplitude spectrum of bathymetric profiles derived 
at various scales by different survey instruments. All data 
were collected on the East Coast continental rise near 
40°27'N, 62°20'W, proximal to the HEBBLE area. Regression 
lines represent the mean bottom roughness model derived from 
fourteen SASS profiles illustrated in Figure 5-13. Model 
predicts the deep-tow derived spectrum almost exactly over 
~1.5 decades of spatial frequency. The spectrum derived 
from a stereo-pair bottom photograph is predicted to within 
0.5 orders of magnitude over ~ decades of spatial frequency 
(200 meters to 2 millimeters wavelength). 
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magnitude is a surprisingly good result. The ability to predict high 
frequency roughness from a low frequency model appears to be quite pos- 
sible, at least for some areas of the world ocean. Ome serious caution 
must be taken into account. The successful prediction from the SASS- 
derived spectrum was only possible by identifying the break in slope at 
3 km. Were data only available to that 3 km wavelength, the lower fre- 
quency estimates would have been used, and prediction error of nearly 
six orders of magnitude incurred in estimating the millimeter scale amp- 
litudes. If it is assumed that such slope breaks are due to changes in 
relief-forming processes at various frequencies, then any estimates of 
high frequency roughness exclusively from a lower frequency model pre- 
sume a continuity of process over all intervening frequencies, and the 


absence of a significant break in the power law form of the spectrum. 
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8. Summary and Conclusions 


A method has been developed to allow a valid stochastic description 
of sea-floor relief to be generated in a. relatively simple statistical 
model. The fundamental statistic used is the amplitude spectrum of spa- 
tial frequency, which is both elementary enough to be generated opera- 
tionally from existing digital bathymetric data bases, and general 
enough to be applied to a variety of engineering applications and scien- 
tific problems. The model allows relatively large areas of the world 
sea floor to be described by as few as two model parameters for simple 
isotropic surfaces. The difficulties of producing a stationary statis- 
tie in a non-stationary environment are minimized with the use of a spa- 
tial domain provineing technique. The model also accounts for the 
directional dependence of anisotropie surfaces. The results of one sim- 
ple experiment indicate that the model may be extrapolated with high 
fidelity to frequencies beyond the resolving capability of surface ship 
sonar systems. This stochastic model when combined with lower frequency 
deterministic models, such as the gridded bathymetric models developed 
by NAVOCEANO, allows a complete description of the sea-floor relief. 

There are numerous avenues available to improve the model. The use 
of ensembles of either data or derived statistics allows the estimation 
errors of the model to be reduced substantially. The availability of 
multibeam sonar systems makes this improvement possible immediately. 
The ensembling of the statistics from the sixteen beams of the SEABEAM 
system would reduce the estimation error to one quarter the results from 
a single beam. Data derived from very narrow beam sonars allow a higher 


rate of sampling, and therefore a wider frequency band, to be available 
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for model generation. The collection of more complete areal bathymetric 
surveys would allow further refinement of the models of sea-floor aniso- 
tropy. Future improvements in deep-towed side-scan sonars and the col- 
lection of additional bottom photographs would allow a quantitative 
estimate to be made of the predictive ability of the model. 

The utility of the model falls into two broad categories; engineer- 
ing application and scientific investigation. The application of the 
model to underwater acoustics is obvious. The model spectrum of the 
surface relief is an important environmental factor in the scattering of 
sound from the sea floor. Efforts are currently underway to develop 
scattering models which utilize such stochastic environmental informa- 
tion. More traditional models require the description of the bottom as 
a faceted surface, which can easily be derived from the derivative spec- 
trum and a knowledge of the probability distribution of depths. Even 
models requiring a fully determined surface can obtain a valid realiza- 
tion by combining the model amplitude spectrum with a randomly generated 
phase spectrum. 

In terms of scientific investigation, the method provides a new 
tool for studying the earth's geological and tectonic processes. The 
resulting patterns of roughness discovered on the Gorda Rise and Oregon 
continental margin illustrate the ability of the method to detect inter- 
esting relationships not obvious by simply studying bathymetric charts. 
Comparing the distribution of roughness on a variety of spreading cen- 
ters, continental margins and other submarine environments, should yield 
valuable insight into the distribution of geologic processes on the sea 
floor. In the case of the patterns on the Gorda Rise, the distribution 


of roughness may well relate to the lava formation processes, which in 
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turn may be related to the distribution of polymetallic sulphide miner- 
als formed in association with hydrothermal vents. 

The ability to quantify the anisotropy of the sea floor also intro- 
duces interesting possibilities for geological investigation. The dis- 
tribution and trend of these measurements provides a new tool for the 
quantitative investigation of deep-sea processes in various environ- 
ments. Also the relationship between sea-floor spreading rate and 
roughness, long assumed qualitatively, could be quantified and analyzed 
in terms of the operative geological processes. 

Finally, the author hopes that this study represents more than sim- 
ply a study of sea-floor roughness. Much effort has been made to pre- 
sent an approach to modelling natural phenomena which could be applied 
to a wide variety of problems in the natural sciences. As most disci- 
plines within the earth sciences represent a marriage of one of the 
“pure” sciences to the study of the earth, this study represents a 
crossing of natural science with engineering statistical methods. Very 
little such work has been done by earth scientists in the past, perhaps 
because the earth is rarely as well ordered as the controlled labora- 
tories of the chemist or physicist. Due to the natural variability of 
the phenomena under study, the geologist’s attempt to describe the earth 
quantitatively is particularly difficult. It is hoped that the approach 
and philosophy presented in this attempt will be of use to future 


investigators. 
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Appendix A.l 


Having concluded that the frequency spectrum of submarine topography 
conforms to a power law functional form, it becomes critical in con- 
structing a model based on this statistic to ensure a valid regression 
fit to the data. Such a regression analysis is necessary both in fit- 
ting the amplitude spectrum itself (amplitude as a function of fre- 
quency) and in fitting the energy envelopes for various spatial fre- 
quency bands for use in spatial-domain provincing. We represent our 
power law function as follows: 

y = f(a,x,b) = ax> where y = dependent variable (amplitude) 

x = independent variable (frequency) 


a,b = regression coefficients 


We can use the property of power law functions that they plot linearly 


on log-log axes and recast the function in log-log space as follows: 
log(y) = log(a) + b log(x) 


which can be easily solved by a simple linear regression. Upon deriving 
the coefficients 6 and log (a), these terms can be reconverted to 


linear-linear space as, 


a 


b=b 


a = antilog (log(a)) 
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to reconstruct our power law form, 


y = ax 


Recalling that the method of least-squares minimizes the total sum 
of squares of residual distance from the observed data to the resultant 
regression curve, the method just described minimizes these distances in 
log-log space. The solution derived in this manner does not minimize 
the total residual distance in linear-linear space, although the esti- 
mates might be quite close. Methods exist for performing the least- 
square fit in either log-log or linear-linear spaces, and these are dis- 
cussed here. The choice of method and the use of weighting schemes 
depend on the distribution of the data being fit as well as the distri- 
bution of the estimation error. In all cases, the error is assumed to 
reside in the amplitude estimates, rather than in the independent 
variable. 

In fitting the energy envelope estimates produced in the delinea- 
tion of stationary provinces (see Chapter 5), the regression must be 
performed in linear-linear space without weighting. The errors asso- 
ciated with the envelope estimates (dependent variable), do not depend 
on frequency band (independent variable) and should not be log-trans- 
formed. In the case of amplitude spectra however, it was show by 
Blackman and Tukey (1958) that estimation error is related to the chi- 
squared distribution and that error bars remain constant in log-log 
space. Under these conditions, the regression analysis is optimally 


performed on log-transformed data. 
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As explained in Scarborough (1930), Article 114., the log-transfor- 
mation of the dependent variable (y) causes the residuals in the least- 
squares residual equations to be of unequal weight. In the case of a 
power law function (which is given as an example in Scarborough (1930) 
and will not be reproduced here), the weighting function is the squared 
dependent variable (y2). Appendix A.2 presents an ASCII FORTRAN 77 pro- 
gram for doing such a weighted regression analysis in log-log space. 
Because of the chi-square distribution of errors associated with the 
amplitude spectrum, however, the residuals in this case are equally 
weighted and no weighting function is required. 

In fitting the envelope estimates of discrete band passes for use 
in “province picking,” there is no requirement for log transformation of 
the data or weighting of residuals. This is the case because the error 
associated with all envelope estimates is theoretically constant. In 
order to derive properly the regression coefficients a and b in linear- 
linear space, it is necessary to use an iterative method. The following 
development is modified from Scarborough (1930), Article 115., to apply 
to the power law functional form. 

We can express the regression coefficients as the sum of initial 


estimates and differences as follows: 


where a,,bo = initial estimates 


4a, db = correction factors 
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It is convenient to use the coefficients derived from performing the 
linear regression on log-transformed data as the initial estimates (a> 
b,) for the iteration process. If we define a new function in terms of 
estimated coefficients as, 


) am by 
y' = £(x,a,,b,) aox 


the discrete values of this approximating function will be, 


y"'y = £(x1,895b 9) 
y'2 = £(x9,895bo) 


¥'n = £(Xq,a9sbo) 


where n is the total number of data pairs used in the regression 


equations. 


Realizing that the coefficients a and b produce the “best fit” solu- 


tions, the minimized residuals are represented as, 


Vo = £(xo,a,b) = Y2 


Ya | £(x, a,b) on Ya 


where Y], Y2»:--ceYq are the observed dependent variables. Substituting 


our approximation yields 
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A = £(x, > ay + Aa, Ls + Ab) Tiyan Dy areveieyy Th 


or 
vs +yi = £(x, » a, + Aa, Le + Ab),i = 1.2,...,n 
Expanding the right side by Taylor's Theorem for the two variables, 


a and b, ytelds 


ae, af, 
v, +9, = £(x,, a, bo) + dalse- » + Ab(sp mp.) > coop $b WpAyooodt 


substituting aay = £(x,5 as? b)) we have 


of, 


af, 
vi tongiriy L Tere + Ab(sp ) tise sitt Lumaladss «ast 


Rearranging terms and dropping higher order derivatives yields 


of, af, 
Vine sags) + Ab) + y' tn ae sb Oo Ap ooonm 


We define 
™,* y's er mw 12 sicleeigD 
where the r's are the residuals for the approximation curve 


y' = £(x 229d 5)- Substituting, we can write our residual equations 
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of, of, 
WA da(s) + Ab(i=—) + 4, 3; 121,2,...,n 
° 


{ ob 1 

Since this system of equations is linear in the correction terms 
4a and Ab, these terms can be derived by the method of least-squares. 

At this point, one variation from the description given by 
Scarborough (1930) is introduced. Because the linearization of the sys- 
tem of equation is only valid for small Aa and Ab, it is possible to 
derive correction terms which yield solutions that extend beyond the 
local neighborhood of linearization. Under these circumstances, it is 
possible that the method will not converge to a valid solution. During 
any particular iteration, this non-convergence would appear as an 
increase in the total sum of squares of the residuals over the previous 
iterations. 

To ensure a decrease in the total residuals (and therefore a con- 
verging solution) during each iteration, the correction terms are scaled 


by a term a to yield, 


= + ad 
a a 


o> @> 


= bb + aAb 


The scaling term a is first set equal to one, and the calculated 
total residuals compared to those calculated in the previous iteration. 
If the residuals do not decrease, the scale factor a is halved until a 
decrease in the total residuals is observed. These new A and b become 
the new “initial estimates” a, and Bo» which are then used in the next 


iteration. The process continues until the values of both Aa and Ab 


reach some suitable minimum, and a final a and b are derived. It is 
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theoretically possible for the solution to converge to a subsidiary min- 
{mum and therefore yield a poor model. However, the selection of the 
initial a, and b, by a regression in log-log space makes convergence on 
subsidiary maxima unlikely. An ASCII FORTRAN 77 subroutine for perform- 
ing this algorithm is presented in Appendix A.3. 

In practice, the regression models for the amplitude spectra (and 
final roughness model) are formulated interactively with a graphic dis- 
play terminal. The operator can interactively edit the bathymetric pro- 
file under examination as well as control the frequency limits included 
in the regression analysis. This allows the software to delete white- 
noise levels, interpolation effects, and other contaminating factors 
from the analysis. The interactive control also allows the operator to 
detect visually spectra composed of two power law segments, and to fit 
each segment individually. This interactive software is included in the 


amplitude spectrum generating software listed in Appendix D. 
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Appendix A.2 


SUBROUTINE POWWGT(A,B,FIRSTX,DELX,Y,N,ILIST,CTOFF1,CTOFF2) 


THIS ROUTINE PERFORMS A BEST FIT TO DATA WITH A POWER LAW 
FUNCTION OF THE FORM Y=A*X**B USING A WEIGHTING METHOD 
AS. DESCRIBED IN SCARBOROUGH(1930), ART. 114. 
INPUTS ARE 
A=COEFFICIENT OF X 
B= EXPONENT OF X 
X= ARRAY OF INDEPENDENT VARIABLE VALUES 
Y= ARRAY OF DEPENDENT VARIABLE VALUES 
N= NUMBER OF DATA PAIRS OF X AND Y 
AMIN= MINIMUM VALUE FOR A CORRECTION TO STOP ITERATION 
BMIN= MINIMUM VALUE FOR B CORRECTION TO STOP ITERATION 


ILIST= 1 FOR SUMMARY OF ITERATION PROCESS, = 0 , NO LISTING 
PROGRAMMED BY C.G.FOX-ADVANCED TECHNOLOGY STAFF ,NAVOCEANO,4/15/83 


DIMENSION X(1024),¥(1024) 


COMPUTE INITIAL ESTIMATE OF A AND B BY PERFORMING A SIMPLE 
LINEAR FIT ON LOG TRANSFORMED DATA 


YSQR=1. 

YSQSUM=0.0 

XPROD=0.0 

XSUM=0.0 

YSUM=0.0 

XSQR=0.0 

X(1)=FIRSTX 

DO 10 I=2,N 

X(1)=X(I-1)+DELX 

WRITE(6,'(80X,2F10.4)') (X(I),¥(1),1=1,20) 
DO 50 J=1,N 

IF(Y¥(J).LE.0.0) GO TO 50 
YTEMP=AL0G10(Y(J)) 

XTEMP=AL0G10(X(J)) 

IF (XTEMP.GT .CTOFF2.0R.XTEMP.LT.CTOFF1) GO TO 50 
YSQR=YTEMP*YTEMP 

XPROD=XPROD+( YTEMP*XTEMP*YSQR ) 
XSUM=XSUM+XTEMP*YSQR 

YSUM=YSUM+YTEMP*YSQR 

YSQSUM=YSQSUM+YSQR 
XSQR=XSQR+(XTEMP*XTEMP*YSQR ) 

CONTINUE 

B=(( YSQSUM*XPROD ) - (XSUM*YSUM) ) /( ( YSQSUM*XSQR ) - (XSUM*XSUM) ) 
A=(YSUM/YSQSUM)-(B*(XSUM/YSQSUM) ) 

A=10.**A 

RETURN 

END 
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Appendix A.3 


SUBROUTINE POWFIT(A,8,X,Y,N,AMIN,BMIN,ILIST) 


THIS ROUTINE PERFORMS A BEST FIT TO DATA WITH A POWER LAW 
FUNCTION OF THE FORM Y=A*X**B USING AN ITERATIVE SUID 
AS DESCRIBED IN SCARBOROUGH(1930), ART. 115. 
INPUTS ARE 
A=COEFFICIENT OF X 
B= EXPONENT OF X 
X= ARRAY OF INDEPENDENT VARIABLE VALUES 
Y= ARRAY OF DEPENDENT VARIABLE VALUES 
N= NUMBER OF DATA PAIRS OF X AND Y 
AMIN= MINIMUM VALUE FOR A CORRECTION TO STOP ITERATION 
BMIN= MINIMUM VALUE FOR B CORRECTION TO STOP ITERATION 
ILIST= 1 FOR SUMMARY OF ITERATION PROCESS, = 0 , NO LISTING 
PROGRAMMED BY C.G.FOX-ADVANCED TECHNOLOGY STAFF ,NAVOCEANO,4/15/83 


DIMENSION X(10),¥(10) 


COMPUTE INITIAL ESTIMATE OF A AND B BY PERFORMING A SIMPLE 
LINEAR FIT ON LOG TRANSFORMED DATA 


XN=FLOAT(N) 
XPROD=0.0 
XSUM=0.0 
YSUM=0. 
XSQR=0. 
RIEL, '(2F10.4)') (X(I),¥(1),121,N) 
DO 50 J=1,N 
IF(Y(J). EQ. 0.0) GO TO 50 
YTEMP=AL0G10(Y(J)) 
XTEMP=AL0G10(X(J)) 
XPROD=XPROD+( YTEMP*XTEMP ) 
XSUM=XSUM+XTEMP 
YSUM=YSUM+YTEMP 
50 XSQR=XSQR+(XTEMP*XTEMP ) 
BO=( (XN*XPROD) -(XSUM*YSUM) ) /( (XN*XSQR ) - (XSUM*XSUM) ) 
AO=(YSUM/XN ) = (BO*(XSUM/XN ) ) 


A0=10.**A0 
COMPUTE SUM OF SQUARES OF THE RESIDUALS 
POLD=0.0 
DO 100 I=1,N 
100 POLD= POLD+((Y(1)-F3(X(I), AO, 80) )**2) 
ITERAT=0 
NBIS=0 
IF (ILIST.EQ.1)WRITE(6,110) 
110 FORMAT(' ITERATION # OF BISECTIONS A B 


*RES IDUALS**2' ) 
IF (ILIST.EQ.1)WRITE(6,120)ITERAT ,NBIS ,AO,80,POLD 
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120 FORMAT(' ',16,10X,14,6X,3(4X,F10.4)) 
ZERO OUT MATRIX TERMS 


150 Al=0.0 
B1=0.0 
D1=0.0 
E1=0.0 
G1=0.0 


COMPUTE TERMS FOR LEAST SQUARES MATRIX CONSTRUCTION 


DO 200 I=1,N 
PARTA=F1(X(I) ,A0,B0) 
PARTB=F2(X(I) ,A0,B0) 
POWF =F 3(X(I) ,A0,B0) 
Al=A1+(PARTA**2) 
B1=B1+(PARTA*PARTB) 
D1=D1+(PARTB**2) 
E1=E1+(PARTA*(Y¥(1I)-POWF )) 
200 Se a eA 
C1=B 


COMPUTE CORRECTION TERMS FOR A AND B 
DIVSOR=(A1*01-B1*C1) 
ACORR=(D1*E 1-B1*G1)/DIVSOR 
BCORR=(A1*GI-C1*E1)/DIVSOR 

CREATE NEW A AND B 


230 A=A0+ACORR 
B=B0+BCORR/AO 


COMPUTE NEW SUM OF SQUARES OF RESIDUALS WITH NEW ESTIMATES 


PNEW=0.0 
00 250 I=1,N 
250 PNEW=PNEW+((Y(1)-F3(X(1I),A,B) )**2) 


TEST FOR CONVERGENT SOLUTION(PNEW < POLD) 
IF NOT, BISECT CORRECTIONS AND RECOMPUTE 


‘EF (PNEW.LT.POLD) GO TO 300 

ACORR=. 5*ACORR 
BCORR=.5*BCORR 

NBIS=NBIS+1 

IF (NBIS.GT.10) GO TO 300 

GO TO 230 


TEST FOR MINIMUM CHANGE OF A AND B 
300 IF(ABS(A-AO).GT.AMIN) GO TO 500 
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IF (ABS (B-B0).GT.BMIN) GO TO 500 
GO TO 900 


C CORRECTION TERM NOT FINE ENOUGH, START NEW ITERATION 


500 ITERAT=ITERAT +1 

A0=A 

BO=B 

POLD=PNEW 

IF (ILIST.EQ.1)WRITE (6,520) ITERAT ,NBIS ,A0,B80,POLD 
520 FORMAT(' ',16,10X,14,6X,3(4X,F10.4)) 

NBIS=0 

GO TO 150 
900 ITERAT=ITERAT +1 

IF (ILIST.EQ.1)WRITE (6,920) ITERAT ,NBIS,A0,B0,POLD 
920 FORMAT(' *,16,10X,14,6X,3(4X,F10.4)) 

RETURN 

END 


FUNCTIONS TO CALCULATE POWER LAW FUNCTION AND PARTIAL 
DERIVATIVES WITH RESPECT TO A AND B 


FUNCTION F1(X3,A3,B83) 

CALCULATE PARTIAL OF A*X**B WITH RESPECT TO A 
F 1=X 3**B3 

A3=A3 

RETURN 

END 


(=) NAOMONOOM 


FUNCTION F2(X4,A4,B4) 
C CALCULATE PARTIAL OF A*X**B WITH RESPECT TO B 
F2=(X4**B4)*ALOG(X4) 
A4=A4 
RETURN 
END 


FUNCTION F3(X5,A5,B5) 
C CALCULATE POWER LAW A*X**B 
F 3=A5*X5**B5 
RETURN 
END 
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Appendix B.1l 


Before generating amplitude spectra from a statistically non-sta- 
tionary sample space, such as the sea floor, one must initially define 
provinces which are relatively homogeneous with respect to the frequency 
spectrum. Since generating an amplitude spectrum directly by Fourier 
transform assumes stationarity over the length of the input series, this 
“province picking” procedure must be performed by estimating the spec- 
trum discretely in the spatial domain. Although the present application 
is new, the concept of estimating spectra in the time/space domain is 
not. Blackman and Tukey (1958) referred to such estimates as “pilot 
spectra” and describe two methods for their calculation. Godfrey (1967) 
also describes a method, very similar to that detailed here, which is 
used for predicting non-stationary time series. This section gives 
details of the procedure used in this study, presents performance tests 
of the algorithm, and include full FORTRAN-77 software for performing 
the analysis. 

The initial step in processing is identical to that required for 
running an FFT. The data must be projected onto a straight-line segment 
and interpolated to even increments in distance. The straight line is 
generated by a simple least-squares fit of the navigation track; the 
best results are obtained when relatively straight line navigation is 
input. Large deviations in the track greatly degrade results. Next, 
points are mapped onto the nearest location on the least-squares line, 
without alteration if they are within a designated “pivot distance,” or 
modified according to the local gradient if beyond this distance. 


Finally, the data are interpolated at a specified interval using a one- 
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dimensional cubic spline, which tends to preserve frequency content. 
The full algorithm is included in SUBROUTINE MAPCTIN. 

The subsequent stage of processing requires band-pass filtering of 
the interpolated data at ten nearly equi-spaced frequency bands. Fil- 
tering is done by sequential application of low-pass and high-pass 
filters, using a non-recursive, symmetric, least-squares filter devel- 
oped by Martin (1957). The frequency response of this bank of filters 
is illustrated in Figure B-l. A review of the performance of the Martin 
filter can be found in McClain and Walden (1979). The filter cutoffs 
are designed to juxtapose at the 100% energy pass level. The total fre- 
quency bank was selected to span .02 - .25 cycles/data interval, or 
wavelengths of approximately 2 - .16 nautical miles for data recorded at 
12 second intervals, or 10 - .8 nautical miles for one minute data. The 
filtering is performed by SUBROUTINE FILTER. 

The next step of the algorithm requires estimating the instanta- 
neous amplitude of all ten band-passed signals. Davis (1974) used a 
simple full-wave rectification, followed by a low-pass smoother to esti- 
mate the energy envelope. For this study, a true Hilbert transform is 
performed and manipulated to generate the energy envelope. The reader 
is referred to Kanasewich (1981) for a full development. Notice that 
the calculation of the envelope via Hilbert transform does require oper- 
ating in the frequency domain. However, because the complete Fourier 
Transform is retained throughout, the requirement of stationarity is not 
applicable. The enveloping algorithm is contained in SUBROUTINE ENVEL. 

The envelope generated in the previous step represents a continuous 
estimate of amplitude through space, for each selected frequency band. 


To estimate the full amplitude spectrum at each point on the profile, 
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% RESPONSE 


FREQUENCY (cycles/data interval) 


Figure B-1 Simplified frequency response of the bank of band-pass fil- 
ters used to estimate amplitude spectra discretely in the 
spatial domain. Actual responses include side lobe leakages 
of less than 5%. 
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all ten amplitude versus frequency estimates are considered. Since the 
power law form of the spectrum is known, this functional form is used in 
an iterative regression procedure described in Appendix A and performed 
in SUBROUTINE POWFIT. The regression coefficients vary widely along the 
profile, and therefore it is desirable, in order to aid in the detection 
procedure, to smooth these estimates. To save calculation time, this is 
done by averaging 91 envelope estimates at each frequency before enter- 
ing the regression routine. This process tends to smear. the boundaries 
somewhat, but makes detection of provinces more reliable. 

With smoothed estimates of the amplitude spectrum available contin- 
uously along a profile, the final stage of processing is to detect sig- 
nificant chances in the estimated spectra and on this basis impose prov- 
inee boundaries. The regression coefficients a and b represent the 
antilog of the y intercept, and the slope, respectively, of the spectrum 
projected in log-log space. The spectral slope, b (which is related to 
the so-called Fractal dimension) is a worthwhile parameter for province 
detection. A simple algorithm is run across the slope estimates to 
detect significant, rapid changes (boundaries). 

The regression coefficient a is correlated to b and therefore does 
mot represent an independent parameter for detection. An alternative is 
to look at the total, band-limited RMS energy of the estimated spectra 
for significant shifts in total energy. These RMS parameters are easily 
calculated using Parseval's Formula (integrating the power spectrum) 
applied to the estimated spectra. In addition to detecting rapid 
changes as with the slope, the RMS is contoured to form segments of some 
minimum size. Although these detectors have proven fairly reliable, it 


is often necessary to interpret certain boundaries where the various 
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detectors disagree. All of the pertinent software is listed in Appendix 
B.2. The program runs with a core size of approximately 540 K bytes, 
and was written in ANSI Standard FORTRAN (1977) for a UNIVAC 1180/2 com- 
puter. The maximum profile length is 4100 interpolated data points and 
145 points are lost from the end of the signal due to filtering. 

Just as an electrical engineer investigates the performance of an 
instrument by operating on signals of known properties, the same tech- 
nique can be used here to test the performance of the province picker 
algorithm. While an engineer might use a step, ramp, or impulse func- 
tion as input, random signals with known spectral forms are used in this 
analysis. Many of the seemingly arbitrary choices of filters, averaging 
procedures, and other design decisions incorporated into the present 
province picker, were selected through feedback from performance tests 
with known signals. 

An obvious choice of a signal with a known amplitude spectrum is 
random “white noise,” with a continuous spectrum of zero slope. There 
are several means of producing such a signal. The simplest is to gener- 
ate pseudo-random noise series of either a normal or uniform distribu- 
tion. Figures B-2 and B-3 illustrate such signals with their spectra, 
which are indeed relatively flat. Notice the large amount of scatter in 
the amplitude estimates. An alternative method of generating “white 
noise” is actually to use a constant amplitude spectrum and uniformly 
distributed random phase spectrum, separate their real and imaginary 
parts, and inverse-transform the signals using the FFI. In this manner 
a perfectly flat, non-varying spectrum is assured, as is illustrated in 
Figure B-4. This is the signal source used in testing for this study, 


and the algorithm is presented in SUBROUTINE WHINOL. 
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Figure B-2 An example of uniformly distributed random noise is shown in 
upper profile. Its corresponding amplitude spectrum (in 
lower profile) shows the flat (zero slope) form indicative 
of white noise. Notice the high degree of scatter in ampli- 
tude estimates. The same computer software was used as that 
used for bathymetric data. 
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Figure B-3 An example of normally distributed random noise in the same 
format as Figure B-2. Notice the tendency of values to 


cluster about the mean (normal distribution) and the large 
scatter of amplitude estimates in the amplitude spectrum. 
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Figure B-4 An example of random “white noise" generated using an 
inverse Fast Fourier Transform. The method of calculation 
forces a perfectly constant amplitude spectrum. This “white 
noise" generator was used for performance testing of the 
province picker algorithm. 
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Figure B-5 illustrates the performance of the province picker when 
white noise is input. Examination of this output results in some inter- 
esting insights into the nature of stochastic processes. Despite the 
known property that the input signal has a perfectly flat, non-variable 
amplitude spectrum, all ten band-pass filtered signals show large fluc- 
tuations in amplitude. The amplitude recorded in the frequency spectrum 
represents simply an average amplitude computed over the length of the 
input time series. At any discrete point, the amplitude at.a frequency 
could be widely removed from the average value. It was discovered 
through experimentation that these fluctuations were damped when wider 
band-pass filters were used; thus, the overlapping filter bank illus- 
trated in Figure B-l was designed. The slope values (plotted above the 
ten band-passed signals) do fluctuate about zero as expected. The 
standard deviation about zero averages about 0.2 over several runs which 
implies (assuming a no.mal distribution) that a change of wsiope of t0.4 
can be detected with 95% confidence. Many of the decisions made in gen- 
erating the regression analysis and smoothing, were designed to minimize 
this fluctuation. 

In order to design a province detector properly, it is necessary to 
combine known signals of differing spectral slopes and RMS energies. 
For the purpose of this study, the resulting signal must also have an 
amplitude spectrum with power law form. One method of generating such a 
signal is through a Markov chain with a probability transition matrix 
which allows subsequent events either to raise or lower a constant 
increment with probability of 0.5 (see Ross, 1980). This signal, which 
is a special case of a random walk, fluctuates about its initial value 


and has a spectral form of 
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Figure B-5 Example of province picker output with "white noise" input 
illustrated in Figure B-4. The input signal is shown at the 
bottom, above which are the output of ten band-pass filters 
convolved with the data. The lowest signal (Band Pass 1) jis 
the lowest frequency pass, while the highest (Band Pass 10) 
is the highest frequency pass. Above each band-passed sig- 
nal is the energy envelope calculated by the Hilbert Trans- 
form method. Notice the large variability in energy for each 
band despite the constant amplitude input. At the top, the 
estimated spectral exponent (slope of log-transformed spec- 
tra) and the band-limited RMS energy calculated along the 
profile are plotted. Standard deviation of the slope par- 
ameter is .2. The algorithm was designed to minimize this 
"natural" variability. 
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It is now possible with two known signals of different spectral 
characteristics to study the response of the province picker as it 
crosses the boundary Secsesn two such signals (provinces). Since it is 
also necessary to detect changes in total RMS energy without a change of 
slope, each signal can be multiplied by a constant using a corollary of 
the addition theorem of Fourier transforms (see Bracewell, 1965). 

Figure B-6 illustrates an artificially generated random signal 
which consists of four distinct provinces. The first half of the signal 
is composed of “white noise" and the second half of 1/s noise generated 
via the random walk technique. Notice the rapid change in the slope 
parameter from 0 to -l at the boundary, which was easily detected. 

Within each half of the signal, the provinces are again divided 
with the second halves (second and fourth quarters) representing a dou- 
bling of the first halves (first and third quarters). Notice the 
obvious change of RMS energy, although the slope parameter is unaf- 
fected, illustrating the importance of detecting on two uncorrelated 
parameters. Several false alarms are observed on the derivative detec- 
tor, but the province boundaries derived from contouring (represented by 
straight lines above the RMS energy profile) correspond to the known 
boundaries in the signal. Notice that provinces one and four show the 
same RMS energy level, but are delineated by their differences in spec- 
tral slope. 

As stated earlier, the automatic detection is an aid to province 
boundary recognition, but in some cases, human intervention is needed to 


resolve inconsistencies. Also, the setting of contour interval or slope 
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Figure B-6 Example of province picker output with mixed input of known 
signals. The first half of the signal is white noise (A = 
as9), joined to the second half of random walk generated 
noise (A = as~') the amplitude of each half is doubled at 
the mid-point. The break in spectral slope parameter is 
easily detected. The change in level of RMS energy is con- 
toured, as shown by straight line segments above the RMS 
energy profile. 
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threshold can =. modified by the investigator depending upon the amount 
of resolution desired. Figure B-7 is a sample profile from real bathy- 
metric data collected by the SASS system in the vicinity of the Gorda 
Rise. Boundaries are less distinct than in the artificially generated 
signal, as one would expect, however distinct changes of RMS energy (and 
in two cases changes of slope) do define quasi-stationary provinces. In 
practice, the slope parameter has offered little independent information 
for province picking, and therefore, construction of a simplified prov- 


ince detector based on RMS energy alone may prove to be adequate. 
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Appendix B.2 


PROGRAM TO SELECT PROVINCES FOR SPECTRUM GENERATION 
PROGRAMMED BY C.G. FOX,ADVANCED TECHNOLOGY STAFF ,NAVOCEANO,5/15/84 


PROVINCE SELECTION IS BASED ON A SPATIAL DOMAIN ESTIMATE OF THE 
FUNCTIONAL DESCRIPTION OF THE SPECTRUM BASED ON A POWER LAW MODEL 


AMPLITUDE= A*FREQUENCY**B 


THE SPECTRUM PARAMETERS A & B ARE ESTIMATED BY PERFORMING A 
REGRESSION ANALYSIS ON THE ENERGY ENVELOPES GENERATED ON TEN 
BAND PASSES OF THE DATA SPACED EVENLY IN FREQUENCY. PROVINCE 
BOUNDARY SELECTION IS BASED ON CHANGES IN THE BAND LIMITED 
RMS LEVEL OF EACH LOCATION ALONG THE TRACKLINE AS DERIVED 

BY INTEGRATING THE RESULTS OF THE REGRESSION MODEL. 


PROGRAM CONTROLS ARE ENTERED IN FREE FORMAT FOLLOWING @XQT 


JTOTL= NUMBER OF TRACKLINES OF DATA TO BE ANALYZED IN THIS ANALYSIS 
MAXIMUM IS CURRENTLY SET TO 50 


IPLOT= 1, PRODUCE PLOT TAPE ON UNIT 10 OF INTERPOLATED DATA, 
OPTIONAL BAND PASS PLOTS, PROVINCES SPECTRAL PARAMETER 
ESTIMATES, AND PROVINCE BOUNDARIES. 

OUTPUT CAN GO TO ZETA OR CALCOMP PLOTTER 
0, NO PLOT TAPE IS PRODUCED 


IBAND= 1, PLOT ALL BAND PASSED SIGNALS WITH ENVELOPES 
0, NO BAND PASSES 
IGNORED IF IPLOT=0 


ISAVE1= 1, SAVE DATA WITHIN EACH PROVINCE ON UNIT 13 FOR LATER 
ANALYSIS IN FFT PROGRAM. 
0, NO DATA SAVED 


ILIST= 1, PRODUCE COMPLETE LISTING OF ANALYSIS 
O, ONLY RUDIMENTARY LISTING IS PRODUCED 


GRID= DESIRED SPACING OF INTERPOLATED DATA IN NAUTICAL MILES 
THIS SPACING IS NORMALLY SLIGHTLY LESS THAN THE ORIGINAL 
SPACING OF THE RAW DATA. 


KPTS= NUMBER OF RAW DATA VALUES TO BE ANALYSED IN EACH PROFILE. 
READING OF THE DATA WILL CEASE WHEN THIS NUMBER IS REACHED 


DATA INPUT IS CURRENTLY DESIGNED SUCH THAT FOLLOWING THIS CONTROL 
CARD, THE PROGRAM EXPECTS A CARD IMAGE CONTAINING THE NAME(LESS 
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THAN 20 CHARACTERS) OF A FILE IN WHICH IS CONTAINED A LIST OF 
ALL FILE NAMES TO BE USED IN THE ANALYSIS. AFTER READING THIS 
FILE NAME, JTOTL FILE NAMES ARE READ FROM THE FILE AND THE 
PROGRAM AUTOMATICALLY OPENS THESE DATA FILES AS NEEDED. USERS 
MUST INSURE THAT THE FILE NAME LIST FILE AND ALL DATA FILES 
ARE AVAILABLE TO THE RUN AT EXECUTION. ONE CONVENIENT WAY 

TO PROVIDE DATA TO THE RUN IS BY CREATING DATA AND LIST 
ELEMENTS IN A PROGRAM FILE AND THEN CREATE AN @ADD ELEMENT TO 
COPY ALL ELEMENTS INTO TEMPORARY FILES BEFORE EXECUTION 


PARAMETER (ISIZ=7000) 
DIMENSION X(ISIZ),Y(ISIZ),Z(ISIZ) ,AX(ISIZ) ,AY(ISIZ),AZ(ISIZ), 
1WGT (240) ,K1(1000) ,K2(1000) ,K3(1000) 
2,FRLOW(10) ,FRHIGH(10) ,ENVEST(10) ,TZ(1SIZ,12) 
3,FRMEAN( 10) 
4, ENERGY (ISIZ),TZ11(1S1Z),1Z12(1S1Z) 
EQUIVALENCE (TZ(1,11),1Z11(1)),(7Z(1,12) ,1Z12(1)) 
CHARACTER*20 FILE(50) 
CHARACTER*20 INFILE 
FRHIGH=HIGH FREQUENCY FOR DATA BANDPASS 
FRLOW= LOW FREQUENCY FOR DATA BANDPASS 
DATA FRLOW/.002, .023,.044, .065, .086,.107,.128, .149,.170,.191/ 
DATA FRHIGH/.041, .062, .083,.104,.125, .146, .167, .188, .209, .230/ 
SET FILTER PARAMETERS FOR BANDPASSES 
FILSLP=.008 
NUMF IL=50 
NFIT IS LENGTH OF AVERAGE OVER ENERGY ENVELOPES 
NFIT=91 
C2MAX=0.0 
SET PARAMETERS TO CONTROL DERIVATIVE PICKS ON SPECTRAL PARAMETERS 
SLWIN=.4 
SLWID1=1. 
SLWID2=14. 
ENWIN=400. 
ENWID1=20. 
ENWID2=100. 
PDEL DETERMINES THE INTERVAL OF RMS ENERGY USED FOR PROVINCE 
DETECTION 
PDEL=.03 
Oat otc NUMBER OF POINTS ALLOWED WITHIN A PROVINCE 
MIN=9 
OSCALE=SCALING FACTOR(INCHES/UNIT) OF ORIGINAL DATA FOR PLOTTING 
OSCALE=-1. 
ESCALE= SCALING FACTOR FOR ENVELOPE UNITS 
ESCALE=.0001 
READ(5,*) JTOTL,IPLOT, IBAND,ISAVE1,ILIST,GRID,KPTS 
IKONT =0 
LTPS=1000 
GRIDKM=GRID*1.852 
IF(IPLOT.NE.1) GO TO 24 
INITIALIZE PLOTTER 
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USE FOLLOWING CALL FOR CALCOMP PLOTTER OUTPUT 
CALL PLOTS(0,0,10) | 
USE FOLLOWING CALL FOR ZETA PLOTTER OUTPUT 
CALL PLOTS(53,0,-10) 
CALL FACTOR(.5) 
CALL PLOT(1.0,1.0,-3) 
READ NAME OF FILE CONTAINING LIST OF DATA FILE NAMES 
24 READ(5,'(A20)') INFILE 
25 OPEN(UNIT=18,FILE=INFILE,STATUS='OLD' ) 
READ LIST OF DATA FILE NAMES FROM FILE INFILE 
00 5 IK=1,JTOTL 
READ(18,'(A20)') FILE(IK) 
5 CONTINUE 
CLOSE(UNIT=18) 
READ X,Y,Z, VALUES FOR THIS TRACK.FIRST POINT IS 
ASSUMED TO BE ORIGIN,PUT END OF FILE AT THE END OF EACH TRACK. 
440 IKONT=IKONT+1 
ZERO OUT TEMPORARY ARRAYS 
DO 441 J=1,10 
DO 441 beens 
441 TZ(1, ge 


442 X(I)= 0. 0 
IAVPLT=0 
INPUT X(I)=LONGITUDE(DEC. DEG.),Y(1)=LAT,Z(1)=DEPTH 
N=# OF POINTS - 1 
IKNT=IKONT 
JPTS=KPTS 
IF(ILIST.EQ.1) WRITE(6,45) FILE(IKONT) 
45 FORMAT(' OPENING FILE ',A20) 
ROUTINE TO READ LAT,LON,DEPTH FROM FILE # IKNT 
CALL PROVRD( JPTS,Y,X,Z,IKNT,FILE) 
4 N=JPTS-1 
IF(ILIST.EQ.1) WRITE(6,'(24H NUMBER OF INPUT POINTS ,I5)")N 
INTERPOLATE DATA ON STRAIGHT LINE 
CALL MAPCTN(GRID,X,Y,Z,N,AX,AY,AZ, JCT) 
eS TE Dee en NUMBER OF INTERPOLATED POINTS ,16)') 
Ci 
839 C=JCT 
AZAVE=0.0 
PLOT INTERPOLATED DATA 
DO 261 I=1,JCT 
261 AZAVE=AZAVE+(AZ(1)/C) 
IF(IPLOT.NE.1) GO TO 26 
XL=(C*.02)+0.5 
DV=50.*GRID 
CALL NEWPEN(1) 


CALL AXES(0.,0.,12HDISTANCE(NM),-12,XL,0.,1.,0.,DV,-1) 
EV=AZAVE-(3.0/0SCALE ) 
DV=1.0/0SCALE 
CALL AXES(0.,-3.,13HORIGINAL DATA,!#,¢.,().,!.,EV,DV,1) 
CALL SYMBOL(XL+.75,-.25,.5,12HINPUT SIGNAL,0.0,12) 
IF(C2MAX.LT.XL) C2MAX=XL 
CALL PLOT(0.0,(AZ(1)-AZAVE)*OSCALE, 3) 
00 26 [=1,JCT 
AJ=I-1 
C=AJ*.02 
D=(AZ(1I)-AZAVE )*OSCALE 
CALL PLOT(C,D,2) 

26 CONTINUE 
00 27 I=1,JCT 

27 Z(1)=AZ(1) 
IF(IPLOT.EQ.1)CALL PLOT(0.0,1.0,-3) 
DO 237 IRUN=1,10 
XRUN=IRUN 
ENVMAX=0.0 
IF(IPLOT.EQ.1.AND.IBAND.EQ.1)CALL PLOT(0.0,4.0,-3) 

C BAND PASS FILTER DATA AT SPECIFIED INTERVAL 
CALL FILTER(AZ,JCT,FRLOW( IRUN) ,FILSLP,NUMFIL,1,WGT) 
CALL FILTER(AZ,JCT ,FRHIGH( IRUN) ,FILSLP ,NUMFIL,O,WGT ) 
IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,3) 
WGTMN=0.0 
XNORM=( (FRHIGH( IRUN)+.2)-FRLOW( IRUN) )*(1./GRIDKM) 
C XNORM=1. 

NTWICE=2*NUMF IL 
DO 31 I=NTWICE,JCT-NTWICE 


Geeel 

IF(IPLOT.NE.1.0R.IBAND.NE.1) GO TO 31 
Gasol 

AJ=1-1 

C=AJ*.02 


B=(AZ(1)/XNORM)*ESCALE . 
IF(I.EQ.NTWICE) CALL PLOT(C,B,3) 
CALL PLOT(C,B,2) 
31 AZ(1)=(AZ(1))/KNORM 
(eS 
IF (IPLOT.NE.1.0R.IBAND.NE.1) GO TO 32 
kasi 
EV=-1.5/ESCALE 
DV=1./ESCALE 
CALL PLOT(0.0,0.0,3) 
GAUIPAXES(OMMONM LHI SL XUN OM SIERO ORR) 
CALL AXES(XL,-1.5,13HFILTERED DATA, -13,3.,90.,1.5,£V,DV,0) 
CALL SYMBOL(XL+.75,-.25,.5,9HBAND PASS,0.0,9) 
CALL NUMBER(XL+5.75,-.25,.5,XRUN, 0.0, -1) 
CALL PLOT(0.0,0.0,3) 
CALL NEWPEN(2) 
C COMPUTE APPROXIMATION TO ENVELOPE OF BAND PASSED DATA 
32 CALL ENVEL(AZ,JCT) 
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44 
836 
837 
840 


841 


845 
846 


850 


IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,3) 
ENVAVG=0.0 


IF(IPLOT.NE.1.0R.IBAND.NE.1) GO TO 321 


DO 44 I=NTWICE,JCT-NTWICE 
AJ=I-1 
C=AJ*.02 
B=AZ(1)*ESCALE 
IF(I.EQ.NTWICE) CALL PLOT(C,B,3) 
CALL PLOT(C,B,2) 
GO TO 322 
DO 44 I=NTWICE,JCT-NTWICE 
AZ(1)=ABS(AZ(I)) 
IF(AZ(1).GT.ENVMAX) ENVMAX=AZ(T) 
ENVAVG=ENVAVG+(AZ(1)/( JCT-2*NUMF IL) ) 
ENVELOPE ENERGY IS STORED IN THE APPROPRIATE COLUMNS 
OF 2-DIMENSIONAL ARRAY TZ 
TZ(1, IRUN)=AZ(1) 
FRDIFF =FRHIGH( IRUN) -FRLOW( IRUN) 
IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,3) 
IF(IPLOT.EQ.1) CALL NEWPEN(1) 
DO 836 [=1,JCT 
AZ(1)=Z(T) 
CONTINUE 
LINEAR REGRESSION OF TEN ENVELOPE ESTIMATES FOR EACH POSITION 
FIRST COMPUTE AVERAGE ENVELOPE 
DO 841 I=1,10 
FRMEAN( I) =(( (FRLOW(I)+.017)+FRHIGH(I))/2.) 
CONVERT CYCLES/DI TO CYCLES/KM 
FRMEAN(I)=FRMEAN(I)*(1./GRIOKM) 
FR1=(FRLOW(1)+.017)*(1./GRIDKM) 
FR2=FRHIGH(10)*(1./GRIDKM) 
DO 850 I=NTWICE , JCT-(NTWICE+NF IT) 
DO 846 J=1,10 
ENVEST(J)=0.0 
FIT=FLOAT(NFIT) 
DO 845 K=1,NFIT 
ENVEST(J)=ENVEST(J)+TZ(I+K, J) 
CONT INUE 
IF(FIT.EQ.0.0) GO TO 846 
ENVEST(J)=ENVEST(J)/FIT 
CONTINUE 
PERFORM ITERATIVE REGRESSION,. STORE B & A IN 
COLUMNS 11 & 12 OF ARRAY TZ 
CALL POWFIT(TZ12(I+((NFIT/2)+1)) ,TZ11(1+((NFIT/2)+1)) SFRMEAN( 1), 
*ENVEST(1),10,0.000001, .01,0) 
PLOT SLOPE PARAMETER B 
IF(IPLOT.NE.1) GO TO 859 
CALL PLOT(0.0,4.0,-3) 
CALL ‘AXES(0.0,0.0,1H ,1,XL,0.,1.,0.,0.,-2) 
CALL NEWPEN(2) 
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CALL AXES(XL,-2.0,5HSLOPE,-5,2.0,90.,1.0,-4.0,2.0,-1) 
CALL PLOT(0.0,0.0,3) 
859 AVSLP=0.0 
KKNT=0 
SLWID=SLWID2-SLWID1 
ENWID=ENWID2-ENWID1 
INWID1=IF IX (ENWID1) 
LSWIDL=IFIX(SLWID1) 
INWID2=IF IX (ENWID2) 
LSWID2=IF IX (SLWID2) 
DO 877 K=(NTWICE+1)+(NFIT/2), JCT-((NTWICE+1)+(NFIT/2)) 
KKNT=KKNT+1 
C=(K-1)*.02 
B=(TZ(K,11))/2. 
AVSLP=AVSLP+(TZ(K,11)/(JCT-(4*NUMF IL+(NFIT+1)))) 
IF(IPLOT.NE.1) GO TO 876 
IF (KKNT.EQ.1) CALL PLOT (C,B,3) 
CALL PLOT(C,B,2) 
876 IF(KKNT.LE.LSWID2) GO TO 877 
IF (K.GE. (JCT-(NTWICE+NFIT/2+1+LSWID2))) GO TO 877 
PREAV=0.0 
POSTAV=0.0 
C DO 8761 KNB=LSWID1,LSWID2 
C ——- PREAV=PREAV+(TZ(K-KNB,11)/SLWID) 
C8761 POSTAV=POSTAV+(TZ(K+KNB,11)/SLWID) 
IF (ABS(POSTAV-PREAV).LT.SLWIN) GO TO 877 
CALL SYMBOL(C,-2.,.2,16,0.0,-1) 
CALL PLOT(C,B,3) 
877 CONTINUE 
C PLOT ARRAY OF RMS ENERGY 
IF(IPLOT.EQ.1) CALL NEWPEN(1) 
KKNT=0 
AVENER=0.0 
DO 878 L=(NTWICE+1)+(NFIT/2) , JCT-((NTWICE+1)+(NFIT/2)) 
CALCULATE BAND LIMITED RMS-FIRST SQUARE FUNCTION AND INTEGRATE 
INTEGRAL ( (A*F **B )**2)=(A**2/(2*B+1)*F **(2*B+1) 
XINTLO=(TZ(L,12)**2/(2.*TZ(L,11)+1.))*FR1**(2.*TZ(L,11)+1. ) 
XINTHI=(TZ(L,12)**2/(2.*TZ(L,11)+1.))*FR2**(2.*TZ(L, 11) +1. ) 
SINCE SPECTRUM IS SYMMETRIC ABOUT ZERO, CAN EVALUATE INTEGRAL 
BY MULTIPLYING EVALUATED RESULT BY TWO 
ENERGY(L)=2.*(XINTHI-XINTLO) 
TO CALCULATE MEAN SQUARE, DIVIDE BY WIDTH 
ENERGY (L)=ENERGY(L)/(2.*(FR2-FR1) ) 
NOW TAKE SQUARE ROOT TO DETERMINE RMS 
ENERGY (L)=SORT (ENERGY (L)) 
878 AVENER=AVENER+(ENERGY (L)/( JCT-(4*NUMFIL+(NFIT+1)))) 
ENSCAL=4.*AVENER ; 
IF(IPLOT.NE.1) GO TO 8782 
eee AXES(0.0,0.0,10HRMS ENERGY,10,2.,90.,0.5,0.0,AVENER/2., 
ee 
CALL PLOT(0.0,0.0,3) 
8782 DO 8781 L=(NTWICE+1)+(NFIT/2) , JCT-((NTWICE+1)+(NFIT/2)) 
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KKNT=KKNT#1 
C=(L-1)*.02 
B=ENERGY(L)*(1./AVENER) 
IF(L.EQ.(NTWICE+1)+(NFIT/2).AND.IPLOT.EQ.1) CALL PLOT(C,B,3) 
IF(IPLOT.EQ.1) CALL PLOT(C,B,2) 
IF(KKNT.LE.INWID2) GO TO 8781 
IF(L.GE.( JCT-( (NTWICE+NFIT/2+1)+INWID2))) GO TO 8781 
PREAV=0.0 
POSTAV=0.0 
DO 8771 KNB=INWID1,INWID2 
PREAV=PREAV+(ENERGY (L-KNB) /ENWID) 
8771 POSTAV=POSTAV+( ENERGY (L+KNB ) /ENWID) 
IF (ABS(POSTAV-PREAV).LT.ENWIN) GO TO 8781 
CALL SYMBOL(C,2.,.2,16,180.0,-1) 
CALL PLOT(C,B,3) 
8781 CONTINUE 
DO 879 I=(NTWICE+1)+(NFIT/2) , JCT-( (NTWICE+1)+(NFIT/2) ) 
RMSSLP=RMSSLP+ABS ((TZ(1I,11)-AVSLP)**2) /( JCT-(4*NUMFIL+(NFIT-1))) 
879 RMSENG=RMSENG+ABS ( (ENERGY (I )-AVENER)**2)/( JCT-(4*NUMFIL+(NFIT-1) )) 
RMSSLP=SQRT(RMSSLP ) 
RMSENG=SQRT (RMSENG ) 
C WRITE(6,'(4(F8.3,2X))') AVSLP,RMSSLP,AVENER,RMSENG 
IST=(NTWICE+1 )+(NFIT/2+1) 
IEND=JCT-(IST-1) 
333 CONTINUE 
C PICK PROVINCE BOUNDARIES ,STORE PROVINCE NUMBER FOR EACH POINT IN X 
DO 34 I=IST,IEND 
X(I)=IFIX(SQRT( ENERGY (1) )/PDEL )+1 
34 CONTINUE 
C OUTPUT POSITIONS Or PROVINCE BOUNDARIES 
C WRITE(06,60) CUT 
C 60 FORMAT(' PROV.NO TYPE 1ST LAT LONG LAST LAT LONG 
C 1 RMS-FILTERED CUT=',F7.5) 
IPNUM=X (IST) 
SUM=0.0 
JST=IST 
KST=1 
DO 35 I=IST,IEND 
JPNUM=X (I) 
SUM=SUM+X (I) 
IC=I-JST 
IF (JPNUM.EQ.IPNUM.AND.I.LT.IEND) GO TO 35 
C IF((IC+1).LT.MIN) GO TO 35 
JJ=I 
C COMPUTE AVERAGE PROVINCE NO. 
AC=IC+1 
AVE=SUM/AC 
IAVE=AVE 
BAVE=IAVE 
AVE=AVE-BAVE 
IF(AVE.GE.0.5) IAVE=IAVE+1 
IPNUM=IAVE 
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K1(KST)=JST 
K2(KST)=JJ 
K3(KST)=IPNUM 
JST=JJ 
SUM=0.0 
IPNUM=X ( JST ) 
IF(KST.EQ.1) GO TO 116 
MST=KST-1 
IF(IAVE.NE.K3(MST)) GO TO 116 
K2(MST )=K2(KST) 
GO TO 35 
116 KST=KST+1 
IF(KST.GT.LTPS) GO TO 117 
35 CONTINUE 
GO TO 118 
117 WRITE(06,119) 
119 FORMAT(' YOU HAVE TOO MANY PROVINCES-INCREASE DIMENSION OR PDEL') 
GO TO 500 
C NOW OUTPUT FINAL PROVINCE BOUNDARIES OF AT LEAST MIN. SIZE 
118 K3(KST)=99 
KST=KST-1 
C SMOOTHE SMALL PROVINCE SEGMENTS 
eae PRVFIX(K1,K2,K3,MIN,KST) 
J2= 
DO 36 [=1,KST 
II=I+1 
IF(K3(II).EQ.K3(I)) GO TO 36 
JST=K1(J2) 
JJ=K2(1) 
IPNUM=K3(J2) 
ISDLAT=AY( JST) 
A=ISDLAT — 
SELAT=ABS (AY ( JST) -A)*60. 
ISOLNG=AX (JST) 
A=ISDLNG 
SELNG=ABS (AX ( JST)-A)*60. 
IDLAT=AY (JJ) 
A=IDLAT 
EELAT=ABS(AY(JJ)-A)*60. 
IDLNG=AX (JJ) 
A=IDLNG 
EMLNG=ABS (AX( JJ)-A)*60. 
IDIFF=K2(1)-K1(1) 
C COMPUTE RMS LEVEL OF HIGH PASSED DATA FOR THIS PROVINCE 
C ALSO AVERAGE SPECTRAL SLOPE AND Y-INTERCEPT 
C3=ABS( JJ-JST+1) 
RMS=0.0 
SLOPEX=0.0 
XINTER=0.0 
DO 37 IE=JST,JJ 
SLOPEX=SLOPEX+(TZ(IE,11)/C3) 
XINTER=XINTER+(TZ(IE,12)/C3) 
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37 RMS=RMS+(ENERGY(IE)/C3) 
C ADJUST ENERGY LEVEL IN LOG-LOG SPACE AFTER NORMALIZATION 
XINTER=10**(ALOG10(XINTER)-.5) 
IPNUM=IF IX(SQRT(RMS )/PDEL )+1 
IF(ILIST.EQ.1)WRITE(6,51)1, 1PNUM, ISDLAT,SELAT, ISDLNG,SELNG, IDLAT, 
1EELAT, IDLNG,EMLNG,RMS , IDIFF , SLOPEX, XINTER 
C LOOP TO INSERT DATA INTO FILE FOR LATER FFT RUNS 
IF(ISAVE1.NE.1) GO TO 49 
ENOMRK=9.999999 
NPTS=JJ-JST 
WRITE(13,46) FILE(IKONT),I ,SLOPEX,XINTER,GRID,NPTS 
WRITE (13,50)1, IPNUM, ISDLAT, SELAT, ISDLNG, SELNG, IDLAT ,EELAT, 
LIDLNG, EMLNG , RMS 
WRITE(13, 47)(AZ(1Q) ,1Q=9ST, Jd) 
XRMS= SRMS (AZ( JST) , NPTS ) 
46 FORMAT(A20,12,1X, F8. 351X693 a XeF 74), 15) 
47 FORMAT ( 10F8. 6) 
WRITE(13,47) ENDMRK 
49 IF (IPLOT.NE.1) GO TO 248 
APNUM=.5+IPNUM*.5 
AJ=JST-1 
BJ=JJ-1 
C1=AJ*.02 
C2=BJ*.02 
CALL PLOT(C1,APNUM, 3) 
CALL PLOT(C2,APNUM, 2) 
50 FORMAT(I5,18,4(1X,14,F6.2),F10.4) 
51 FORMAT(15,18,4(1X,14,F6.2),F10.4,16,2F10.6) 
248 J2sII 
36 CONTINUE 
880 IF(IKONT.EQ.JTOTL) GO TO 500 
C REPOSITION PLOTTER FOR NEXT PLOT 
IF (MOD(IKONT,4).NE.0) GO TO 882 
C2MAX=C2MAX+15 .0 
IF(IPLOT.EQ.1) CALL PLOT(C2MAX,-44.0, -3) 
C2MAX=0.0 
GO TO 440 
882 IF(IPLOT.EQ.1) CALL PLOT(0.0,8.0,-3) 
GO TO 440 
500 IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,999) 
IF(ISAVE1.GT.0) CLOSE(UNIT= 13) 
STOP 
END 


157 


READER SUBROUTINE FOR PROVPICK 
SUBROUTINE PROVRD(NP,XLAT ,XLONG, DEPTH, IKONT FILE) 
DIMENSION XLAT (3000) ,XLONG( 3000) , DEPTH ( 3000) 
CHARACTER*20 FILEN 
CHARACTER*20 FILE(20) 
KNTER=0 
400 FILEN=FILE(IKONT ) 
WRITE (6,433 )FILEN 
433 FORMAT(' OPENING FILE ',A20) 
OPEN(UNIT=12,FILE=FILEN,STATUS='OLD' ,FORM='FORMATTED' ) 
DO 420 I=1,NP 
438 READ(12,434,END=470) XLAT(1I),XLONG(I),IDEPTH 
434 FORMAT (9X ,F12.8,F13.8,17) 
KNTER=I 
IF(IDEPTH.LE.0) PRINT *, "DATA POINT SKIPPED' 
IF (IDEPTH.LE.0) GO TO 438 
DEPTH(I )=IDEPTH*. 0018288 
420 IF(I.LT.30)WRITE(6,*)XLAT(I) ,XLONG(I) ,DEPTH(T) 
GO TO 490 
470 NP=KNTER 
490 CLOSE (UNIT=12) 
WRITE (6,492 )NP 
492 FORMAT(' NUMBER OF POINTS READ =',15) 
RETURN 
END 
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SUBROUTINE MAPCTN(GRID,X,Y,Z,N,AX,AY,AZ, ICT) 
ROUTINE TO MAP DATA VALUES(Z) WITH ASSOCIATED POSITIONS 
(X=LONG DECIMAL DEG.,Y=LAT DECIMAL DEG. )FROM A RELATIVELY STRAIGHT 
SEGMENT OF SURVEY TRACK ONTO A STRAIGHT LINE,ADJUST THE (Z) VALUE FOR 
THE AMOUNT OF SHIFT REQUIRED BY THE MAPPING AND INTERPOLATE NEW (Z) 
VALUES (AZ) AT AN EQUALLY SPACED DISTANCE(GRID IN DECIMAL NAUTICAL MI. 
AND ASSOCIATED POSITIONS AX,AY ALONG THIS STRAIGHT LINE. 
C** N=NO.OF ORIGINAL X,Y,Z INPUT PTS., ICT= NO.OF OUTPUT PTS. 
C** NOTE-SINCE THIS IS A CARTESIAN MAP WHICH DOES NOT ACCOUNT FOR LONG 
C CONVERGENCE ,IT SHOULD NOT BE USED FOR SEGMENTS COVERING A LARGE 
C LATITUDE RANGE. 
C** NOTE-ORIGINAL XYZ DATA IS DESTROYED IN THIS ROUTINE 
C*** THIS ROUTINE CALLS GINT, SPLINE ,SPLICO,SORTY 
DIMENSION X(7000),Y (7000) ,Z( 7000) ,AX (7000) ,AY( 7000) ,AZ( 7000), 
*DIST (7000) 
00 55 I=1,7000 
55 DIST(I)=0.0 
JOIR=1 
IF (ABS (X(N)-X(1)).GE.ABS(Y(N)-¥(1))) JDIR=0 
ATER=9999.99 
RAD=0. 00029089 
BLONG=X (1)*60.0 
BLAT=Y(1)*60.0 
DO 1 I=1,N 
IF (JDIR.EQ.1)GO0 TO 3 
X(1)=BLONG-X (1)*60.0 
Y(1)=¥(1)*60.0-BLAT 
GO TO 1 
3 AY(1I)=BLONG-X (1)*60.0 
X(1)=¥(1)*60.0-BLAT 
Y(1I)=AY(T) 
1 CONTINUE 
AN=N 
C MAP INPUT POSITIONS ONTO STRAIGHT LINE TO PREPARE DATA FOR 
C INTERPOLATION ON AN EQUAL DISTANCE BASIS 
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5 D=0+¥(1)*xX(1) 

Al=(C*B-D*A) / (AN*B-A**2) 

A2=(C-A1*AN)/A — 

IF (ABS (A2).LT.0. 000001 )A2=0. 000001 

LEAST SQUARES LINE IS Y=A1+A2X 

NOW SEARCH FOR A POINT LESS THAN PIVOT DISTANCE FROM TRACK LINE 

TO USE FOR 1ST PIVOT AND MAP PTS.ONTO TRACK. WITH CORRECT Z VALUE 
C POINTS-LT PIVOT DIST FROM TRACK WILL HAVE THEIR POSITS MAPPED ONTO 
C We aees BUT THEIR Z VALUE WILL NOT BE CHANGED 

UM=0.0 
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DO 6 I=1,N 
DIST(1)=ABS((¥(1)-A2*X(1)-A1)/ (SQRT(A2**2+1.0))) 
6 SUM=SUM+DIST(1)**2 
PIVOT =SQRT (SUM/AN ) 
C REMOVE POINTS THAT MAY HAVE A BAD POSITION 
XP IVOT=3. O*PIVOT 
J=0 
DO 42 I=1,N 
IF (DIST(1).GT.XPIVOT) GO TO 42 
J=J+1 
X(J)=X(T) 
Y(J)=¥(1) 
Z(J)=Z(1) 
DIST(J)=DIST(I) 
42 CONTINUE 
N=J 
AN=N 
00 7 I=1,N 
IF (DIST(I)- PIVOT) 8,8,7 
7 CONTINUE 
8 A3= -1.0/A2 
AX (1 )=(A1+A3*X (1 )-¥(1))/(A3-A2) 
AY(I)= A2*AX(1) +Al 
AZ(1)=Z(1) 
IA=I+1 
C NOW WORK BACKWARDS ON TRACK TO PICK UP POINTS THAT FAILED 
C PIVOT TEST 
11 J=I-1 
LEO) 1254259 
9 DELZ=(Z(J)-Z(1))/SQRT((X(1)-X (J) )**2 +(¥(1)-¥ (0) )**2) 
AX (J )=(A1+A3*X (J )-¥(d)) /(A3-A2) 
AY (J )=A2*AX (J)+A1 
AZ (J )=(DELZ*SQRT ( (AX (1 )-AX (J) )**2+ (AY (1)-AY(J))**2))+AZ(T) 
I=J 
GO TO ll 
C NOW WORK FORWARD ON TRACK TO PICK UP REMAINING PTS. 
12 DO 13 I=IA,N 
IF (DIST(I)- PIVOT) 14,14,15 
14 AX(1)=(A1+A3*X (1 )-¥(1))/(A3-A2) 
AY(T)= A2*AX(I) +Al 
AZ(1)=Z(T) 
GO TO 13 
15 J= I-1 
DELZ=(Z(1)-Z(9))/SQRT((X(1)-X(J))**2 +(¥(1)-¥(0))**2) 
AX (1 )=(A1+A3*X (1)-¥(1))/(A3-A2) 
AY(1)=A2*AX(I) +Al 
AZ (1 )=(DELZ*SQRT ( (AX (1 )-AX (J) )**2+ (AY (1 )-AY (J) )**2) )+AZ (0) 
13 CONTINUE 
SORT INPUT SO THAT THE INDEPENDENT VARIABLE IS MONOTONIC AND REMOVE 
CLOSELY SPACED POINTS TO CONTROL SPLINE INTERPOLATION 
DAVE =0.0 
DO 44 [=2,N 
J=I-1 
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44 DAVE=DAVE+X (1)-X (J) 
DAVE =0. 3* (ABS (DAVE /AN ) ) 
CALL SORTY(AY,AX,AZ,Y,X,Z,N,1,1,DAVE ) 
IF(JDIR.LT.1) GO TO 92 
DO 93 I=1,N 
Z(1)=AZ(T) 
X(1)=AY(1) 
93 Y(1)=AX(1) 


0 
92 DO 95 


94 CONTINUE 
RESET ORIGIN TO FIRST MAPPED POINT 

BLAT =BLAT +Y¥ (1) 

BLONG=BLONG-X (1) 

DX=X (1) 

DY=Y(1) 

DO 71 I=1,N 

X(1)=X (1 )-DX 
71 Y(1)=¥(1)-DY 

NOW INTERPOLATE ALONG TRACK AT DESIRED GRID SPACING 
COMPUTE APPROX.LENGTH OF LONGITUDE FOR THIS TRACK 
USE APPROX. MIDLATITUDE AS BASIS AND TABLE 6(BOWDITCH) 

NN=N/2 

AL =(Y(NN )+BLAT )*RAD 

A=(111415.13*COS (AL )-94.55*COS (3. *AL )+.012*COS(5.*AL ))/60.0 

UNITS OF A ARE METERS/MINUTE OF LONGITUDE 
COMPUTE APPROX.NAUTICAL MILES/MINUTE OF LONGITUDE 
19 B=A*(1.0/1852.0) 
CONVERT X COORDINATE OF MAPPED POSITION TO MILES AND COMPUTE DISTANCE 
DOWN TRACK. 

00 25 I=1,N 
25 DIST(1)=SQRT ((X(1)*B)**2+Y (1 )**2) 

IF (ABS (X(N)).LT.0.00001) X(N)=0.00001 

THETA=ATAN2(Y(N), (X(N)*B) ) 

CALL GINT(GRID,N,0.0,DIST(N), ICT, DIST,Z,AX,AZ) 
STORE NEW POSIT OF INTERPOLATED PT.IN AX(LONG,AY(LAT) DECIMAL DEG. 
NEW INTERPOLATED VALUE OF Z IS STORED IN AZ 

DO 26 I=1,ICT 

AY (1T)=(AX (T)*SIN (THETA )+BLAT )/60.0 
26 AX (I )=(BLONG-(AX (I )*COS(THETA)/B) )/60.0 

RETURN 

END 
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C FUNCTION TO CONVERT LATITUDE INTO MERIDIONAL PARTS: 
FUNCTION YMP(Z) 
DATA AP/0.7853981634/ 
Y=ABS (Z)*0. 290888209E -03 
T=TAN (AP+Y*0. 5) 
YM=7915. 7045*ALOG10(T )-23. 268932*SIN(Y) 
YMP=YM*SIGN(1.0,Z) 
RETURN 
END 


162 


SUBROUTINE SORTY(X,Y,Z,AX,AY,AZ,K,KODE , JCODE ,GRID) 


C Y=INPUT VARIABLE TO BE SORTED, X,Z=VALUES ASSOCIATED WITH Y 
C K=LENGTH OF Y,IF KODE=1,VALUES OF Y WHICH ARE WITHIN 1 GRID INT OF 
C PREVIOUS VALUE ARE REMOVED,IF KODE=0 ALL VALUES OF Y ARE 
C RETAINED, IF JCODE=+1,Y IS SORTED IN INCREASING ORDER,IF JCODE=-1, 
C Y IS SORTED IN DECREASING ORDER,OUTPUT IS SORTED VALUES OF 
C Y WITH ASSOCIATED X AND Z 
DIMENSION Y¥(20),X(20),Z(20),AX (20) ,AY(20) ,AZ(20) 
KB = K 
CODE =JCODE 
Jzl 
1295 alsa el 
JCT=0 
AY(J) = Y(T) 


132 TEMP= CODE*(AY(J)-Y(I+1)) 
IF ((ABS (TEMP) )-GRID+0.0001) 122,120,120 
TZO SIR GUEME) lal 36; 123 
2 eles a te 
Wao 6 dh)) Ss) obeys ite y2e ac) 
123) eh =) 2 
AY(J) = 
AX(J) = 
AZ(J) = 
RU a ee 


122 IF(KODE) 120,120,136 
136 KD=1+2 
IF (KD-KB) 124,124,139 
139 K=K-1 
KB=KB-1 
GO TO 125 
124 DO 126 JD=KD,KB 
JF = J0-1 
Y(JD) 
X(JD) 
Z(J0) 
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125 IF(JCT) 127,127,128 
Weil Nid) YG) 


AX(J) = X(1) 
aa) 2 OU) 
KT = 2 


128) dander 
IF(J - K) 131,133,133 
131 DO 134 KA = KT,KB 
JT = KA - 
Y(JT) 
X(JT) 
IY. Aaa) 
KB = KB - 1 
GO To 129 
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133 IF(JCT) 137,137,138 
137 KB=KB+1 
138 AY(K) 


135 


AX (K) 
AZ(K) 
DO 135 
Y(T) 
X(1) 
Z(1) 
RETUR 
END 
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Y(KB- 1) 
X(KB- 1) 
Z(KB- 1) 
I = 1,K 
AY(1) 
AX(1) 
AZ(1) 
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SUBROUTINE GINT(DELX,M,XBGN,XEND,ICT,X,Y,AX,AY) 
C MODIFICATION OF ORIGIONAL GINT SO THAT INPUT DATA IS NOT DESTROYED 
C GENERAL 1-D SPLINE INTERPOLATION FOR MIN.STORAGE 
C INTERPOLATES OVER 50 INPUT PTS WITH OVERLAP 
C DELX=DESIRED INTERPOLATION INTERVAL,X AND Y ARE INPUT WITH X=INDEP. 
C VARIABLE,AX AND AY ARE OUTPUT,M=LENGTH OF X, XBGN AND XEND=DESIRED 
C BEGINNING AND ENDING VALUES OF AX,ICT=LENGTH OF AX 
C***NOTE IF M MOD 50 IS LESS THAN 2 INTERPOLATED OUTPUT MAY STOP SHORT 
C OF XEND 
C THIS ROUTINE REQUIRES SPLINE AND SPLICO SUBROUTINES 
DIMENSION X(1),Y(1),AX (1) ,AY(1) 
ICT=(ABS (XEND-XBGN ) /DELX )+1.0 
ISECT=M/50 
IA=ISECT*50 
IC=0 
JJ=0 
MM=49 
IF(IA.EQ.0) MM=M 
IF((M-IA).GE.2) JJ=1 
C OVERLAP INTERPOLATION INTERVALS BY 2 INPUT PTS SO SPLINE ROUTINE 
C IS DIMENSIONED TO MAX OF 53 PTS 
JCONT=1 
KA=1 
KC=MM+1 
IF (IA.EQ.0) KC=KC-1 
K=1 
53 IKT= (ABS (X(MM)-XBGN)/DELX )+1.0 
IF (MM.EQ.M) IKT=ICT 
ATER=9999.999 
DO 26 I=KA,IKT 
AJ=1-1 
X INT =Ad*DELX+X BGN 
IF (XINT.GT.XEND) GO TO 58 
CALL SPLINE (X(K),Y(K) ,KC,XINT, YINT ATER) 
AX (1) =X INT 
26 AY(I)=YINT 
JCONT =JCONT +1 
IF(IA.EQ.0) GO TO 56 
IF(JCONT.GT.ISECT) GO TO 55 
JC=JCONT*50 
57 IMM=MM-1 
K=IMM 
MM=MM+50 
IF(IC.EQ.1) MM=M 
KA=IKT+1 
KC=JC-IMM+1 
GO TO 53 
55 IF(JJ.LT.1) GO TO 56 
IC=1 
JJ=0 
JC=M 
GO TO 57 
58 ICT=I-1 
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RETURN 

56 ICT=IKT 
RETURN 
END 
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SUBROUTINE SPLINE (X,Y,M,XINT,YINT,ATER) 
SEE PENNINGTON REF. FOR DESCRIPTION OF THIS SUBROUTINE 
DIMENSION X(1),Y(1),C(4,53) 
K=0 
IF (X(1)+Y¥ (M)+Y¥ (M-1)+X (M-1)+Y (M-2)-ATER) 10,3,10 
10 CALL SPLICO (X,Y,M,C) 
ATER= X(1)+Y¥(M)+Y(M-1)+X (M-1)+Y (M-2) 
K=1 
3 IF(ABS(XINT-X(1)).LT.0.00001) GO TO 1 
IF (XINT-X(1)) 70,1,2 
70 K=1 
GO TO 7 
1 YINT=¥(1) 
RETURN 
2 IF (ABS(XINT-X(K+1)).LT.0.00001) GO TO 4 
IF (XINT-X(K+1))6,4,5 
4 YINT=Y(K+1) 
RETURN 
5 K=K+1 
IF(M-K) 71,71,3 
71 K=M-1 
GO TO 7 
6 IF (ABS(XINT-X(K)).LT.0.00001) GO TO 12 
IF (XINT-X(K))13,12,11 
12 YINT=Y(K) 
RETURN 
13 K=K-1 
GO TO 6 
11 YINT=(X(K+1)-XINT )*(C(1,K)*(X(K+1)-X INT )**2+C(3,K) ) 
YINT=YINT+(XINT-X (K) )*(C(2,K)*(XINT-X (K) )**2+C(4,K)) 
RETURN 
7 PRINT 101,XINT 
7 CONTINUE 
101 FORMAT(' CAUTION VALUE AT POSITION',F10.2,'WAS EXTRAPOLATED' ) 
GO TO ll 
RETURN 
END 
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SUBROUTINE SPLICO (X,Y,M,C) 
DIMENSION X(1),¥(1),C(4,53) ,D(53) ,P(53) ,E{53),A(53,3),B(53) ,Z(53) 

MM=M-1 

DO 2 K=1,MM 

D(K) =X (K+1)-X (K) 

P(K)=D(K)/6. 
E(K)=(Y(K+1)-¥(K))/D(K) 

00 3 K=2,MM 

B(K)=E(K)-E(K-1) 
A(1,2)=-1.-D(1)/D(2) 
A(1,3)=0(1)/D(2) 
A(2,3)=P(2)-P(1)*A(1,3) 
A(2,2)=2.*(P(1)+P(2))-P(1)*A(1,2) 
A(2,3)=A(2,3)/A(2,2) 
B(2)=B(2)/A(2, 2) 

DO 4 K=3,MM 

A(K,2)=2.* (P(K-1)+P(K) )-P(K-1)*A(K-1, 3) 
B(K)=B(K)-P(K-1)*B(K-1) 
A(K,3)=P(K) /A(K, 2) 

B(K)=B(K) /A(K, 2) 

Q=D(M-2)/D(M-1) 
A(M,1)=1.+0+A(M-2, 3) 
A(M,2)=-Q-A(M, 1)*A(M-1, 3) 
B(M)=B(M-2)-A(M, 1)*B(M-1) 
Z(M)=B(M)/A(M, 2) 

MN =M-2 

DO 6 I=1,MN 

K=M-I 

Z(K)=B(K)-A(K, 3)*Z(K+1) 
Z(1)=-A(1,2)*Z(2)-A(1,3)*Z(3) 

DO 7 K=1,MM 

Q=1./(6.*D(K)) 

C(1,K)=Z(K)*Q 

C(2,K)=Z(K+1)*Q 

C(3,K)=¥(K) /D(K)-Z(K)*P(K) 

C(4,K)=¥(K+1)/D(K)-Z(K+1)*P(K) 
RETURN 

END 
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SUBROUTINE AXES(X,Y,IBCD,NC,AXLEN,ANG,DELTIC ,FIRSTV,DELVAL ,NDEC) 


C MODIFIED CALCOMP AXIS SUBROUTINE ---RANKIN,NOV.1971 
C X,Y COORDINATES OF STARTING POINT OF AXIS IN INCHES 
C IBCD AXIS TITLE 
C NC NUMBER OF CHARACTERS IN TITLE 
C IF NC IS(-),HEADING AND TICKS ARE BELOW THE AXIS 
C AXLEN FLOATING POINT AXIS LENGTH IN INCHES 
C ANG ANGLE OF AXIS FROM HORIZONTAL IN DEGREES 
C DELTIC DISTANCE BETWEEN TIC MARKS IN INCHES 
C FIRSTV SCALE VALUE AT FIRST TIC MARK 
C DEL VAL SCALE INCREMENT 
C NDEC NUMBER OF DECIMAL PLACES OF TIC ANNOTATION PLOTTED(PUNCH 
C -1 IF ONLY INTEGER(NO DECIMAL POINT)IS DESIRED) 
C -2 IF NO ANNOTATION DESIRED IE. ONLY TICKS 

DIMENSION IBCD(10) 

A=1.0 

KN=NC 

IF(NC)1,2,2 

1 A=-A 
KN=-NC 


2 XVAL=FIRSTV 
STH=ANG*0. 0174533 
CTH=COS (STH ) 
STH=SIN(STH) 
OXB=-0.1 
DYB=0.15*A-0.05 
XN=X+DXB*CTH-DYB*STH 
YN=Y+DYB*CTH+DXB*STH 
NT IC=AXLEN/DELTIC+1.0 
NT=NTIC/2 
00 10 I=1,NTIC 
IF(NDEC.EQ.-2) GO TO 12 
C CHANGED NUMBER HEIGHT FROM .105 TO .15 
CALL NUMBER(XN,YN,0.15,XVAL ,ANG,NDEC ) 
12 XVAL=XVAL+DELVAL 
XN=XN+CTH*DELTIC 
YN=YN+STH*DELTIC 
IF (NT )10,11,10 
11 Z=KN 
OXB=-0.07*Z+AXLEN*0.5 
DYB=0. 325*A-0.075 
XT =X+DXB*CTH-DYB*STH 
YT =Y+DYB*CTH+DXB*STH 
C CHANGED HEIGHT FROM .14 TO .18 
CALL SYMBOL (XT,YT,0.18,1BCD(1) ,ANG,KN) 
10 NT=NT-1 
CALL PLOT (X+AXLEN*CTH, Y+AXLEN*STH, 3) 
IF(NDEC.EQ.-2) GO TO 14 
OXB=-0.07*A*STH 
DYB=0.07*A*CTH 
GO TO 13 
14 DXB=-0.05*A*STH 
DYB=0.05*A*CTH 
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13 


20 


10 
20 


A=NTIC-1 

XN=X +A*CTH*DELTIC 
YN=Y+A*STH*DELTIC 

00 20 [=1,NTIC 

CALL PLOT(XN, YN, 2) 

CALL PLOT(XN+DXB, YN+DYB, 2) 
CALL PLOT(XN,YN, 2) 
XN=XN-CTH*DELTIC 
YN=YN-STH*DELTIC 

CONT INUE 

RETURN 

END 

FUNCTION SRMS(X,N) 
DIMENSION X(1) 

AVE=0.0 

SRMS=0.0 

DO 10 I=1,N 

AVE =AVE +(X(1)/FLOAT(N)) 
00 20 I=1,N 

SRMS =SRMS + ((X(1)-AVE )**2) 
SRMS=SRMS /FLOAT (N ) 

SRMS =SQRT (SRMS ) 

RETURN 

END 
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SUBROUTINE FILTER(X,NP,CUT,H,N,K,WGT ) 
C GENERAL ROUTINE TO HIGH OR LOW PASS A SET OF EQUALLY 
C SPACED DATA USING MARTIN FILTERS. 
C X=INPUT DATA AND OUTPUT DATA,NP=NO.OF POINTS IN X,CUT=NORMALIZED 
C CUTOFF OF FILTER IN CYCLES/DATA INTERVAL ,H=SLOPE PARAMETER, 
CN=HALF LENGTH OF FILTER,TOTAL LENGTH=2N+1, 
C K=1=HIGH PASS,=0 FOR LOW PASS. 
C WEIGHTS STORED IN WGT 
DIMENSION X(1),WGT(1) 
NA=N+1 
WGT (1)=2.0* (CUT+H) 
C CENTER WEIGHT STORED IN LOCATION 1 
SUM=0.0 
PI=3. 1415926 
DO 61 I=2,NA 
P=I-1 
Q=1.0-16. O*H**2*P**2 
IF (ABS(Q).GT.0.0001) GO TO 62 
WGT (I )=SIN(2.*PI*P*(CUT+H))/(4. O*P) 
GO TO 61 
62 WGT(1)=((COS(2.*PI*P*H))/Q)*((SIN(2.*PI*P*(CUT+H)) )/(PI*P)) 
61 SUM=SUM+WGT (I) 
DELTA=1.-(WGT(1)+2.*SUM) 
AX=2*N4+1 
IF(K.LT.1) GO TO 78 
DO 65 I=2,NA 
65 WGT(1)=(WGT (1)+DELTA/AX )*(-1.0) 
WGT (1)=1.0-(WGT (1)+DELTA/AX ) 
GO TO 79 
78 DO 80 I=1,NA 
80 WGT(I)=WGT(I)+DELTA/AX 
79 NB=NP-N 
C CONVOLVE WEIGHTS WITH DATA. 
DO 63 I=NA,NB 
II =I+1-NA 
SUM=0.0 
DO 64 J=1,NA 
Jl=I+J-1 
J2=I-J+1 
64 SUM=SUM+WGT (J)*(X(J1)+X (J2)) 
63 X(II)=SUM-WGT (1)*X (I) 
C SHIFT FILTERED DATA TO CORRECT LOCATION AND ZERO ENDS. 
II =NB+1-NA 
DO 67 I=1,I1 
J=II+1-1 


67 X(L)=X( 
I 
68 X(I)=0. 
DO 69 I=NC,NP 

69 X(I)=0.0 


RETURN 
END 
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SUBROUTINE ENVEL (DATA,NP) 
THIS ROUTINE RECEIVES A TIME SERIES OF LENGTH NP - 
IN ARRAY DATA, AND RETURNS ENVELOPE ESTIMATES BASED 
ON A HILBERT TRANSFORM METHOD AS DESCRIBED IN KANASEWICH 
(1981) ,P. 362-368. 
C.G.FOX ,ADVANCED TECHNOLOGY STAFF ,NAVOCEANO-5/12/83 
DIMENSION DATA(1) 
COMPLEX XDATA(2048) 
NSTART=1 
NPSEG=NP 
IF (NPSEG.LE.0) RETURN 
PERFORM ENVELOPE CALCULATIONS IN SEGMENTS OF 2048 POINTS 
IF (NPSEG.GT.2048) THEN 
NPSEG=NPSEG-2048 
LENGTH=2048 
GO TO 100 
ELSE 
LENGTH=NPSEG 
NPSEG=0 
END IF 
TRANSFER INPUT DATA TO REAL PORTION OF COMPLEX ARRAY XDATA 
00 110 I=1,LENGTH 
XDATA(I )=(1.0,0.0)*DATA(NSTART+(I-1)) 
ZERO FILL ARRAY IF LESS THAN 2048 POINTS 
IF (LENGTH.EQ. 2048) GO TO 200 
DO 120 I=LENGTH+1,2048 
XDATA(T )=(0.0,0.0) 
FOURIER TRANSFORM TO FREQUENCY DOMAIN USING FFT 
CALL NLOGN(11,XDATA,-1.0) 
COMPUTE HILBERT TRANSFORM IN THE FREQUENCY DOMAIN. THE 
TRANSFORM IS COMPUTED AS I*SGN(OMEGA)*F (OMEGA), WHERE 
OME GA=FREQUENCY F(OMEGA)=FOURIER TRANSFORM 
SGN (X )=1, (X>0), =0, (X=0) , =-1(X<0) 
THIS OPERATION INDUCES A NINETY DEGREE PHASE SHIFT ON 
COMPLEX PLANE. THE CALCULATION IS DONE BY ZEROING X(1), 
(OMEGA=0) ,TRANSFERRING REAL TO IMAGINARY,MINUS IMAG TO REAL 
FOR X(2)-X(N/2),(OMEGA>O), AND TRANSFERRING MINUS REAL TO 
IMAGINARY,AND IMAGINARY TO REAL FOR X(N/2+1) TO X(N), 
(OMEGA<O ) 
XDATA(1)=(0.0,0.0) 
DO 250 [=2,1024 
TEMPR=REAL (XDATA(I)) 
TEMP I =AIMAG(XDATA(1)) 
XDATA(I )=CMPLX (-TEMPI, TEMPR ) 
DO 260 [=1025,2048 
TEMPR=REAL (XDATA(T )) 
TEMPI=AIMAG(XDATA(I)) 
XDATA(I )=CMPLX (TEMPI , -TEMPR ) 
INVERSE TRANSFORM SHIFTED DATA 
CALL NLOGN(11,XDATA,+1.0) 
CALCULATE ENVELOPE AFTER EQUATION 21.3-1 OF KANASEWICH 
DUE TO SYMMETRIES IN THE TRANSFORM, IMAGINARY PORTION 
OF XDATA IS NEAR ZERO AND CAN BE IGNORED 
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DO 300 I=1,2048 
N=NSTART+(I-1) 
300 _DATA(N)=SQRT ( (DATA(N)*DATA(N ) )+(REAL (XDATA(I))**2)) 
NSTART=NSTART +2048 
GO TO 10 
END 
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SUBROUTINE NLOGN(N,X,SIGN) 
C*****NLOGN COMPUTES THE DISCRETE FOURIER TRANSFORM BY THE FAST FOURIER 
C*****TRANSFORM METHOD.REFERENCE ROBINSON,PAGE 63. 
C*****xTHIS SUBROUTINE CONTAINS 51 STATEMENTS. 
C*****N=POSITIVE INTEGER FOR THE POWER OF 2 DESIRED. 
C*****X=INPUT AND OUTPUT DATA ARRAY (COMPLEX ). 
C*****STGN=-1.0 FOR FOURIER TRANSFORM,+1.0 FOR INVERSE FOURIER TRANSFORM. 
DIMENSION M(15),X(2) 
COMPLEX X,WK,HOLD,Q 


EX =2e=N 
DO 11 =1,N 

TMC) = 2** (NT) 
00 4L = 1,N 
NBLOCK = 2**(L -1) 
LBLOCK = LX / NBLOCK 
LBHALF = LBLOCK / 2 
K =0 
DO 4 IBLOCK = 1,NBLOCK 
FK =K 
FLX = LX 


V = SIGN * 6.2831853071796 * FK / FLX 
WK = CMPLX(COS(V),SIN(V)) 

ISTART = LBLOCK * (IBLOCK - 1) 

00 2 I = 1,LBHALF 


J = ISTART + I 
JH = J + LBHALF 
Q = X(JH) * WK 


X(JH) = X(J) - Q 
X(J) = X(J) + Q 


2 CONTINUE 
00 3 I = 2,N 
II =I 
IF(K.LT.M(I)) GO TO 4 
3K =K - MI) 
4K =K + MII) 
K = 0 
DON Jes a leniex 
IF (K.LT.J) GO TO 5 
HOLD = X(J) 
X(J) = X(K + 1) 
X(K + 1) = HOLD 
5 00 6 I =1,N 
II =I 
IF(K.LT.M(I)) GO TO 7 
6 K =K - M(I) 
7X =K + M(II) 
IF (SIGN.LT.0.0) RETURN 
DO 8 I = 1,LX 
8 X(I) = X(I) / FLX 
RETURN 
END 
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SUBROUTINE PRVFIX(K1,K2,K3,MIN,KST ) 
DIMENSION K1(1000) ,K2(1000) ,K3(1000) ,11(200) ,T2(200) ,13(200) 
SUBROUTINE FOR USE WITH THE PROVINCE PICKER PROGRAM 
SEARCHES PROVINCE BOUNDARIES FOR SEGMENTS LESS THAN MIN 
SEGMENTS ARE INCORPORATED INTO ADJACENT PROVINCES IF THEY 
ARE EQUAL-SETS BOUNDARY TO INCLUDE SEGMENT IN ROUGHEST 
PROVINCE IF ADJACENT PROVINCES ARE UNEQUAL-DELETES 
SEGMENTS AT ENDS OF TRACKLINE 
C.G.FOX-LDGO-JULY 6,1982 


MINSAV=MIN 
MIN=0 
10 ISKIP=0 
ICYCLE=0 
MIN=MIN+5 
IF (MIN.GT .MINSAV )MIN=MINSAV 
DO 100 I=1,KST 
IE GLEYCRESGE=)G0) 10/195 
IX=I-ISKIP 
Il =I+l 
TEST FOR ADJACENT PROVINCES OF SAME NUMBER 
IF(K3(1).NE.K3(II)) GO TO 30 
T1(IX )=K1(1) 
T2(IX )=K2(1T) 
T3(IX )=K3(1) 
ICYCLE=1 
ISKIP=ISKIP+1 
GO TO 100 
TEST FOR MINIMUM PROVINCE SIZE 
30 IDIFF=K2(I )-K1(1) 
IF(IDIFF.LT.MIN)GO TO 40 
NORMAL PROVINCE - STORE IN T ARRAY AND CONTINUE 
T1(IX )=K1(1) 
T2(IX )=K2(1) 
T3(IX )=K3(1) 
GO TO 100 
LESS THAN MINIMUM SIZE-TEST NUMBER OF ADJACENT PROVINCES 
IF FIRST OR LAST PROVINCE OF LINE, DELETE 
40 IF(I.NE.KST.OR.IX.NE.1) GO TO 50 
ISKIP=ISKIP+1 
GO TO 100 
TEST IF SURROUNDING PROVINCES ARE EQUAL 
50 IF(T3(1X-1).NE.K3(II)) GO TO 70 
IF EQUAL, EXTEND SECOND BOUNDARY OF PREVIOUS PROVINCE 
TO END OF FOLLOWING PROVINCE 
T2(1X-1)=K2(11) 
ICYCLE=1 
ISKIP=ISKIP+2 
GO TO 100 
ADJACENT PROVINCES ARE DIFFERENT-SET BOUNDARY CLOSER TO 
LOWER VALUED PROVINCE 
70 IF(T3(1X-1).LT.K3(11)) GO TO 75 
T1(1X )=K1(T) 
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T2(1X )=K2(IT) 

T3(1X )=K3(1) 

IF ((K2(11)-K1(11)).GT.IDIFF )T3(1X )=K3(11) 

GO TO 80 

T2(1X-1)=K2(1) 

T1(IX )=K1 (IT) 

T2(IX )=K2(11) 

T3(1X )=K3(1) 

IF ((K2(1I1)-K1(11)).GT.IDIFF )T3(1X )=K3(11) 
BISECT SEGMENT AND PUT HALVES INTO ADJACENT PROVINCES 
K=(K2(1)+K1(1))/2 
T2(1X-1)=K 
T1(IX )=K 
T2(1X)=K2(11) 

T3(IX )=K3(1T) 
ISKIP=ISKIP+1 
ICYCLE=1 

GO TO 100 
ICYCLE=ICYCLE -1 
CONTINUE 

IF (ISKIP.EQ.0.AND.MIN.EQ.MINSAV) GO TO 500 
KST=KST-ISKIP 
DO 150 I=1,KST 
K1(I)=T1(1) 
K2(T)=T2(1) 
K3(I )=T3(1) 
K3(KST +1 )=99 

G0 TO 10 

RETURN 

END 
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FUNCTION DECDEG(I,X ) 
C CONVERT DEGREES+MINUTES TO DECIMAL DEGREES 
DECDEG=ABS(I)+X/60. 
DECDEG=SIGN(DECDEG,1) 
RETURN 
END 
SUBROUTINE DEGMIN(A,J,B) 
C CONVERT DECIMAL DEGREES(A) TO DEGREES(J)+MINUTES(B) 
J=INT(A) 
B=ABS ( (A-FLOAT (J))*60. ) 
RETURN 
END 


177 


oy a 4 digi’ - 


ah i fh } 


age 


a.t ee sei pase 


“s ENV CE mi 
ath) TAM ne (AYE 


Lae TEL fon nan) ns 


Appendix C.1 

The addition of two sinusoids of different amplitude and phase, but 
identical wavelength results in another sinusoid of the same wavelength, 
but amplitude and phase which is a function of the amplitudes and phases 
of the component sinusoids. Assuming two sinusoids of differing ampli- 
tudes (vy and Vp) and phases (84 and 9,), their sum can be expressed in 
terms of a new sinusoid of amplitude v, and phase 96 as 

vi cos(6-8,) + vps cos(8-8,) = ie ° cos(9-8,) (Cc-1) 


Using the addition formula of sinusoids, 


VN cos9 ° cos8 , Tvia. sind ° sind, + VR’ cosd ° coso, + 


VR’ sinO « sino, = Vone cos9 ° cos8, TaVon sind ° sin, (C=-2) 


Combining terms, and equating sine and cosine terms 


(Jr, cos8 , +v,° cos® 5) * cos = (vo . cos6 ,) * cosd (C=3) 


and 


(vy ° sind, + vg ° sin?) ¢ sind = (ve ° sin®¢) ¢ sind (c=-4) 


Since C-3 and C-4 are true for any value of 9, the unknown v, and 9, 


¢c 


are determined from: 
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v, ° cos6, + vp ° cos®g = vg * cos8> (c-5) 
Va ° sind, + vg ° sind, = Vo° sin8> (C-6) 


Squaring both sides and summing equations C-5 and C-6 removes 


the 9, terms (cos26 + sin*0 = 1); 


23 vac ° cos”6 , + vR- ° cos “0, +2e Via 2 Vp? cos6 , e cos6, + 


2 


MC 


2 e 2 e 2 e e e e - 
sin*9, + vp sin“9, + 2° v, * vg * sind, * sind, (C-7) 


VN A 


Again using the relationship cos26 + sin6 = 1, the addition 
formula for cosines, and taking the root of the result gives the 


formula for the Vc term as, 
1 
Vo = (v4? + vB aa ONAN SANGO (cos(8,-0,))/2 (C=8) 


The unknown resulting phase angle 9¢ can be derived by dividing 


(c-6) by (C-5), 


Vince sind, + VR ° sin, 
tano, = (C-9) 


Vino? sin® , + VR ° cos, 


or 


2 Va ° sin® , + VR ° sin?, 
Q, = tan ( ) (C-10) 


Va ° cos8 , + vp ° cos8, 
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These general results show the dependence of a sinusoid of given 
wavelength on a linear combination of two other sinusoids of the same 
wavelength. The same combination process could be extended to any number 
of component sinusoids. 

Appendix C.2 contains FORTRAN -77 software to perform an iterative 


fit of a sinusoid to data. 
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Appendix C.2 


SUBROUTINE COSFIT(A,B,PHI,X,Y,N,AMIN, BMIN, PHIMIN, ILIST) 


THIS ROUTINE PERFORMS A BEST FIT TO DATA WITH A TRIGONOMETRIC 
FUNCTION OF THE FORM Y=A+B*(COS(4*PI*(X-PHI))) USING AN ITERATIVE 
METHOD AS DESCRIBED IN SCARBOROUGH(1930), ART. 115. 
INPUTS ARE 
A=CONSTANT 
B=AMPLITUDE OF SINUSOIDAL COMPONENT 
PHI= ANGLE OF MAXIMUM AMPLITUDE 
X= ARRAY OF DEPENDENT VARIABLE VALUES 
Y= ARRAY OF DEPENDENT VARIABLE VALUES 
N= NUMBER OF DATA PAIRS OF X AND Y 
AMIN= MINIMUM VALUE FOR A CORRECTION TO STOP ITERATION 
BMIN= MINIMUM VALUE FOR B CORRECTION TO STOP ITERATION 
PHIMIN= MINIMUM VALUE FOR PHI CORRECTION TO STOP ITERATION 
ILIST= 1 FOR SUMMARY OF ITERATION PROCESS, = 0 , NO LISTING 
PROGRAMMED BY C.G.FOX-ADVANCED TECHNOLOGY STAFF ,NAVOCEANO, 8/30/83 


DIMENSION X(1024),Y(1024) 
REAL 19,J9,K9,L9 


COMPUTE INITIAL ESTIMATE OF A AND B BY PERFORMING A SIMPLE 
LINEAR FIT ON LOG TRANSFORMED DATA 


XN=FLOAT(N) 
XPRCD=0.0 
XSUM=0.0 
YSUM=0.0 
XSQR=0.0 
WRITE(6,'(2F10.4)') (X(I),¥(1),1=1,N) 
FIND MEAN LEVEL AND CONVERT X TO RADIANS 
DO 50 I=1,N 
X(T)=X(1)/57.2957795 
50 AO=A0+(Y(1)/XN) 
LOCATE LARGEST POSITIVE DIFFERENCE FROM THE MEAN AND ITS PHI 
00 60 I=1,N 
IF((Y(I)-AO).LT.BO) GO TO 60 
BO=Y(1I)-AO 
PHIO=X(1) 
60 CONTINUE 


COMPUTE SUM OF SQUARES OF THE RESIDUALS 


POLD=0.0 

DO 100 I=1,N 
100 POLD=POLD+((Y(1)-F3(X(1),A0,BO,PHIO) )**2) 

ITERAT=0 

NBIS=0 

IF (ILIST.EQ.1)WRITE (6,110) : 
110 FORMAT(' ITERATION # OF BISECTIONS A B 
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*PHIO RES IDUALS**2) 
IF(ILIST.EQ.1)WRITE (6,120) ITERAT,NBIS,A0,BO,PHIO,POLD 


120 FORMAT(’ ',16,10X,14,6X,4(4X,F10.4)) 


150 A9=0. 


200 


ZERO OUT MATRIX TERMS 


B9=0. 
c9=0. 
E9=0. 
F9=0. 
19=0. 
J9=0. 
K9=0. 
L9=0. 


(efookokoko koe) 


COMPUTE TERMS FOR LEAST SQUARES MATRIX CONSTRUCTION 


DO 200 I=1,N 
PARTA=F1(X(I),A0,B0) 
PARTB=F2(X(1I),A0,PHIO) 
PARTP=F4(X(1),A0,80,PHIO) 
POWF =F 3(X(I),A0,B0,PHIO) 
A9=A9+(PARTA**2 ) 
B9=B9+(PARTA*PARTB ) 
C9=C9+(PARTA*PARTP ) 
E9=E9+(PARTB**2) 
F9=F9+(PARTB*PARTP ) 
19=19+(PARTP**2 ) 

J9=J9+ (PARTA* (Y(1)-POWF ) ) 
K9=K9+ (PARTB*(Y(1I)-POWF ) ) 
L9=L9+(PARTP*(Y(1)-POWF ) ) 
D9=B9 

G9=C9 

H9=F9 


COMPUTE CORRECTION TERMS FOR A AND B 


D=A9*E9*19+B9*F 9*G9+C9*D9*H9 -C9*E 9*G9-B9*D9*19-A9*F9*HO 

ACORR=(( (E9*19-F9*H9 )*J9 )+ ( (C9*H9-BO*I9)*KI)+( (BO*F9-E9*C9)*LO))/D 
BCORR=( ( (G9*F9-D9*I9 )*J9)+( (A9*19-G9*C9)*K9)+ ( (D9*C9-AI*F9)*LI9))/D 
PCORR=(( (D9*H9-G9*E9 )*J9 )+ ( (G9*B9-A9*H9 )*K9)+( (AI*E9-B9*DI)*LO) )/D 


CREATE NEW A AND B 


230 


A=A0+ACORR 

B=B0+BCORR 

PHI=PHIO+PCORR 

COMPUTE NEW SUM OF SQUARES OF RESIDUALS WITH NEW ESTIMATES 


PNEW=0.0 
DO 250 I=1,N 
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250 PNEW=PNEW+((Y(1)-F3(X(I),A,B,PHI) )**2) 


TEST FOR CONVERGENT SOLUTION(PNEW < POLD) 
IF NOT, BISECT CORRECTIONS AND RECOMPUTE 


IF(PNEW.LT.POLD) GO TO 300 
ACORR=.5*ACORR 
BCORR=.5*BCORR 
PCORR=.5*PCORR 

NBIS=NBIS+1 

IF(NBIS.GT.10) GO TO 300 
GO TO 230 


TEST FOR MINIMUM CHANGE OF A AND B 


300 IF(ABS(A-AO).GT.AMIN) GO TO 500 
IF (ABS(B-BO).GT.BMIN) GO TO 500 
IF (ABS(PHI-PHIO).GT.PHIMIN) GO TO 500 
GO TO 900 


CORRECTION TERM NOT FINE ENOUGH, START NEW ITERATION 


500 ITERAT=ITERAT+1 
AO=A 
BO=B 
POLD=PNEW 
IF (ILIST.EQ.1)WRITE (6,520) ITERAT ,NBIS,A0,B80,PHIO,POLD 
520 FORMAT(' ',16,10X,14,6X,4(4X,F10.4), 
NBIS=0 
GO TO 150 
900 ITERAT=ITERAT+1 
IF (ILIST.EQ.1)WRITE (6,920) ITERAT,NBIS,A0,B0,PHIO,POLD 
920 FORMAT(' ',16,10X,14,6X,4(4X,F10.4)) 
RECONVERT RADIANS TO DEGREES 
DO 930 I=1,N 
930 X(1)=X(1)*57.2957795 
PHI=PHI*57.2957795 
PRINT *,A,B,PHI 
RETURN 
END 


FUNCTIONS TO CALCULATE POWER LAW FUNCTION AND PARTIAL 
DERIVATIVES WITH RESPECT TO A AND B 


FUNCTION F1(X3,A3,B3) 
cae PARTIAL OF FUNCTION WITH RESPECT TO A 
Fl=1. 

A3=A3 

X3=X3 

B3=B3 

RETURN 
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END 


FUNCTION F2(X4,A4,PHI4) 
CALCULATE PARTIAL OF FUNCTION WITH RESPECT TO B 
F2=COS(2.*(X4-PHI4) ) 
A4=A4 
RETURN 
END 


FUNCTION F3(X5,A5,B5,PHI5) 
CALCULATE FUNCTION 

F 3=A5+B5*COS (2.* (X5-PHI5) ) 

RETURN 

END 


FUNCTION F4(X6,A6,B6,PHI6) 
CALCULATE PARTIAL WITH RESPECT TO PHI 
F4=2.*SIN(2.*(X6-PHI6) ) 
A6=A6 
B6=B6 
RETURN 
END 
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Appendix D 


NOTE: Many of the subroutines required by the following 
routines are common to the province selection 
software listed in Appendix B.2. 


GENERATE AMPLITUDE SPECTRUM OF DEPTHS OUTPUT FROM PROVINCE PICKER. 
PROGRAM IS A SLIGHT MODIFICATION OF THE ORIGINAL CREATED BY 
T.M.DAVIS, IN WHICH THE INPUTS HAVE BEEN SET WITHIN THE CODE 
AND ONLY THOSE PARAMETERS NECESSARY FOR DAILY USE HAVE BEEN LEFT 
FOR INPUT BY THE USER. INPUT IS EXPECTED FROM UNIT 13, WHICH IS 
THE UNIT NUMBER USED BY PROVPICK. THE PROGRAM ADDS AN INTERACTIVE 
GRAPHICS SECTION WHICH ALLOWS MODIFICATION OF THE INPUT SERIES, 
DELETION OF SECTIONS, AND CONTROL OF THE INDEPENDENT VARIABLE FOR 
REGRESSION ANALYSIS. ROUGHNESS MODEL PARAMETERS A & B ALONG WITH 
LATITUDE AND LONGITUDE INFORMATION FOR EACH PROFILE ARE OUTPUT 
TO UNIT 14 FOR LATER USE IN PROVCHART PROGRAM. 


PROGRAMMMED BY C.G.FOX,ADVANCED TECHNOLOGY STAFF ,NAVOCEANO 5/17/84 
THE FOLLOWING COMMENTS ARE FROM THE ORIGINAL PROGRAM BY T.M.DAVIS 


PROGRAM TO COMPUTE PREWHITENED,CORRECTED AND SMOOTHED AMPLITUDE 
SPECTRUM CUT,H,N=PARAMETERS TO CONTROL PREWHITENING AND SMOOTHING 
FILTERS ,ANORM=SPECTRUM NORMALIZATION IN TERMS OF DATA INTERVALS 
XORG,YORG=LOG-LOG PLOT ORIGIN IN POWERS OF 10, ITG=HEADING 
IUNIT=5 IF CONTROL DATA ON CARDS = TAPE UNIT IF ON TAPE 
SAME FOR JUNIT FOR INPUT DATA 

TAPE UNIT 2=QUTPUT 
CODE IFEOF=1 IF EOF FOLLOWS EACH SET 

FIRST SET OF FILTER PARAMETERS ARE FOR HIGH PASS PREWHITENING 
SECOND SET ARE FOR LOW PASS SMOOTHING 

IF N(1)=-1 COSINE TAPER IS APPLIED,=0 NO TAPER OR PREWHITENING 
LEAVE HIGH OR LOW PASS PARAMETERS BLANK IF NO FILTER DESIRED 
INPUT DATA LIMIT IS 2048 PTS,NP=NO.OF INPUT PTS 

IF ANORM IS BLANK SET TO 1 IN PROGRAM 

USES SUBROTINE NLOGN FROM ROBINSON 

IN FIRST CONTROL CARD NSETS=NO.OF SETS OF DATA. THIS RUN 

SET KPHA=1 IF PHASE SPECTRUM IS DESIRED 

IORGN=NO.OF INPUT DATA PT.TO USE AS ORIGIN FOR PHASE SPECTRUM 

IF PHASE IS DESIRED CODE IPHA=0 IF PHASE IN DEGREES OR CODE IPHA 
=NO.OF DATA INTERVALS/INCH FOR PLOT IF YOU WANT PHASE IN DATA INT. 
SET JCODE=1 IF ONLY PLOT OF FINAL SMOOTHED SPECTRUM DESIRED 

SET JCODE=-1 FOR NO PLOT AT ALL 

CODE ILIST=1 IF YOU DESIRE ONLY LISTING OF SMOOTHED SPECTRUM 


**TMPORTANT-IF YOU HAVE ALREADY ADDED ZEROES TO THE BEGINNING OR END 


OF YOUR INPUT DATA AND YOU REQUEST PREWHITENING OR COSINE TAPER 
YOU MUST CODE IADJ=1 FOR PROPER EXECUTION 

IMEAN=1 IF YOU WANT THIS MEAN REMOVED 
DIMENSION DATA(4500) ,STO(2500) ,CUT(2) ,H(2) ,N(2) ,ADATA( 4100) 
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1,SDATA( 4500) , BDATA( 2050) ,CDATA( 2050) , PHASE (2050) ,SMPH( 2050) 
COMMON /BOUND/I, IPNUM, ISDLAT ,SELAT, ISDLNG,SELNG, IDLAT,EELAT, 
LIDLNG,EMLNG 
CHARACTER*4 ITH(8) ,ITG(8) 
CHARACTER*1 RESPON 
C READ(5,* )NSETS , JCODE , IUNIT, JUNIT, IFEOF ,GRIDNM 
NSETS=1000 
JCODE =1 
IUNIT=5 
JUNIT=2 
IFEOF=1 
C GRID SPACING IS READ FROM DATA FILE 
IF(JCODE.GE.0) CALL PLOTS(STO,2500,2) 
CALL FACTOR(.65) 
JSET=1 
IF(JCODE.GE.0) CALL PLOT(1.0,0.0,-3) 
READ(IUNIT ,*)NP ,KXORG,KYORG,ANORM1, (CUT(I),H(I),N(1),1=1,2),IMEAN, 
1 KPHA,IPHA, IORGN, ILIST, IADJ 
NUMBER OF POINTS IS READ FROM DATA FILE, MAXIMUM IS HELD TO 2000 
NPTS=2000 
KXORG=-2 
KYORG=-7 
ANORM1=2.0 
CUT(1)=-1.0 
H(1)=.2 
N(1)=3 
CUT (2)=.08 
H(2)=.2 
N(2)=3 
IMEAN=1 
KPHA=0 
IPHA=0 
IORGN=1 
ILIST=0 
IADJ=1 
245 CONTINUE 
APHA=IPHA 
XORG=10.0**KXORG 
YORG=10. 0**KYORG 
IOUT =7 
IF (ABS (ANORM1) .LT.1.0E-20) ANORM=1.0 
CALL PRVOUT(NP,DATA, ITG,XINTER,SLOPEX ,GRIDNM) 
GR IOKM=GRIDNM*1. 852 
C XINTER=10** (ALOG1O(XINTER )-(ALOG10(1. /GRIDKM)+.5)) 
2451 CTOFF2=.5*(1./GRIDKM) 
CTOFF 1=.01*(1. /GRIDKM) 
CTOFF 1=AL0G10(CTOFF 1) 
CTOFF 2=AL0G10(CTOFF 2) 
CT ISV=CTOFF 1 
CT 2SV=CTOFF 2 
C SAVE DATA ARRAY IN SDATA ARRAY 
DO 2452 I=1,NP 
2452 SDATA(I)=DATA(I) 


(on ap ap) 
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25 IF (ANORM1.GT.1. )ANORM=1./FLOAT (NP) 
38 IF(CUT(1).NE.-1.0) GO TO 29 
FUNDAMENTAL FREQUENCY ,H(1)=.2,N(1)=3 
CUT (1)=.8/ (FLOAT (NP ) ) 
ACL) =<12 
N(1)=3 
IF (ANORM1.GT.1.) ANORM=1./(FLOAT(NP)-(2*(N(1)-1))) 
29 IF(IADJ.NE.1) GO TO 383 
COMPUTE INDEX OF FIRST NON-ZERO POINT 
28 00 385 J=1,NP 
Kl=J 
IF (DATA(J).NE.0.0) GO TO 386 
385 CONTINUE 
WRITE (IOUT, 82) XORG, YORG,ANORM, ITG 
82 FORMAT(15H PLOT ORIGIN X=,E11.4,4H Y= ,£11.4,15HNORMALIZATION =, 
1£11.4,30A1) 
COMPUTE INDEX OF LAST NON-ZERO POINT 
386 DO 391 J=NP,K1,-1 
K2=J 
IF (DATA(J).NE.0.0) GO TO 392 
391 CONTINUE 
SHIFT DATA TO LEFT AND RECOMPUTE NO.OF INPUT POINTS 
392 NP=K2-K1+1 
DO 393 J=1,NP 
K=K1+J-1 
393 DATA(J)=DATA(K) 
DEMEAN DATA AND APPLY COSINE TAPER 
383 SUM=0.0 
XNP=NP 
DO 388 J=1,NP 
388 SUM=SUM+DATA( J) 
AVE =SUM/XNP 
IF(N(1).GE.0.AND.IMEAN.EQ.0) GO TO 381 
DO 389 J=1,NP 
389 DATA(J)=DATA(J)-AVE 
CALL BATHXS(DATA,NP,GRIDKM) 
IF(N(1).GE.0) GO TO 381 
DO 382 J=1,NP 
AM=J-1 
TNP=NP -1 
TNP=TNP/2.0 
AM=3.14159* (AM-TNP ) /TNP 
382 DATA(J)=DATA(J)*0.5*(1.0+COS (AM) ) 
RAISE NO.OF PTS TO A POWER OF 2(M)AND STORE IN INP 
381 KOUNT=1 
DO 5 M=1,12 
[=2**™ 
IF(I-NP)5,6,7 
5 CONTINUE 
6 INP=I 
GO TO 9 
7 INP=I 
NPP=NP+1 
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DO 8 J=NPP, INP 
8 DATA(J)=0.0 
9 AINP=INP 
DELR=1.0/AINP 
DELK=(1.0/AINP) /GRIDKM 
C INP= TOTAL NO.OF DATA PTS.TO A POWER OF 2 
c DELR=NORMALIZED FREQUENCY INTERVAL IN CYCLES PER DATA INTERVAL 
WRITE (IOUT ,50)DELR, INP,NP 
50 FORMAT(' FREQ. INT=',£8.3,' POWER OF 2 DATA PTS=',16,' NO.OF INPUT 
1 DATA PTS=',16) 
WRITE ( IOUT, 560) IORGN 
560 FORMAT(' ORIGIN FOR PHASE SPECTRUM IS INPUT PT.NO.',14) 
C NOW STORE DATA IN COMPLEX FORM IN ADATA 
53 00 2. J=1,INP 
1=(J*2)-1 
ADATA(I)=DATA(J) 
IA=I+1 
2 ADATA(IA)=0.0 
c COMPUTE AND PLOT NORMALIZED AMPLITUDE SPECTRUM 
CALL NLOGN(M,ADATA, -1.0) 
INX=INP+1 
IF (KOUNT-2) 301,301, 302 
C CORRECT FFT OF PREWHITENING FILTER FOR PHASE SHIFT 
302 ANN=N(1) 
CALL TSHIFT(M,ADATA, ANN, DELR) 
DO 303 I=1,INX,2 
JX=(1+1)/2 
J=I+1 
BDATA( JX )=(SQRT (ADATA(I)**2+ADATA(J)**2) ) *ANORM 
IF(ADATA(I).LT.0.0) BDATA( JX )=-BDATA( JX ) 
303 CONTINUE 
GO TO 204 
301 DO 201 I=1,INX,2 
JX=(1+1)/2 
J=I+1 
201 BDATA( JX )=( SQRT(ADATA(I)**2 +ADATA(J)**2) )*ANORM 
c COMPUTE ROUGH PHASE SPECTRUM FROM PREWHITENED DATA 
G OR COSINE TAPERED DATA 
IF (KPHA.EQ.1.AND.N(1).LE.0) GO TO 567 
IF (KPHA.EQ.1.AND.KOUNT.£Q.2) GO TO 567 
GO TO 568 
567 ANN=IORGN-1 
C CORRECT PHASE SPECTRUM FOR DESIRED ORIGIN 
CALL TSHIFT(M,ADATA,ANN, DELR) 
DO 569 I=1,INX,2 
JX=(I+1)/2 
B=JX-1 
B=B*DELR*360.0 
J=I+1 
IF (ABS (ADATA(I)).LT.1.£-20) ADATA(I)=ADATA(I)+1.€-20 
PHASE ( JX )=(ATAN2(ADATA(J) ,ADATA(I)) )*57.295779 
569 IF(IPHA.GT.0.AND.B.GT.0)PHASE ( JX )=PHASE ( JX )/B 
IF(IPHA.GT.0) GO TO 541 


92 


93 
94 


95 


97 


83 
99 


206 


WRITE (IOUT, 570) 
FORMAT (' ROUGH PHASE SPECTRUM(DEG)' ) 
GO TO 542 
WRITE (IOUT , 543) 
FORMAT(' ROUGH PHASE SPECTRUM(DATA INTERVALS )') 
IF(ILIST.NE.1) WRITE (IOUT, 81) (PHASE (J),J=1,JX) 
IF (KOUNT -2) 202,203, 204 
WRITE (IOUT, 75) 
FORMAT (28H AMP.SPECT.OF ORIGINAL DATA ) 
IF(ILIST.NE.1) WRITE (IOUT,81)(BDATA(J),J=1,Jx ) 
FORMAT (10E11.4) 
GO TO 206 
WRITE (IOUT, 76) 
FORMAT (30H AMP.SPECT.OF PREWHITENED DATA) 
GO TO 77 
WRITE (IOUT,91) CUT(1),H(1),N(1) 
FORMAT(' SPECT.OF PREWHITENING FILTER WITH PHASE REVERSALS CUTOF 
PSY Pe G Se Bs SP Alek ats2 =',13) 
DO 92 J=1,JX 
BDATA(J) =BDATA(J)/ANORM 
IF(ILIST.NE.1) WRITE (IOUT, 81)(BDATA(J),J=1, JX) 
CORRECT SPECTRUM OF FILTERED DATA FOR PREWHITENING FILTER 
D0 93 J=1,JX 
IF (ABS (BDATA(J)).LT.1.£-20)BDATA(J)=BDATA(J)+1.E-20 
BDATA(J)=ABS (CDATA(J)/BDATA(J)) 
WRITE (IOUT , 94) 
FORMAT(' FINAL ROUGH CORRECTED SPECTRUM WITHOUT PHASE REVERSALS') 
IF(ILIST.NE.1) WRITE (IOUT,81)(BDATA(J),J=1, Ux ) 
GO TO 206 
CONT INUE 
SMOOTH FINAL ROUGH CORRECTED SPECTRUM 
COMPUTE LOW PASS WEIGHTS 
K=2 
IF(N(2).EQ.0) GO TO 330 
STORE ROUGH SPECTRUM IN DATA 
00 97 J=1,JX 
DATA(J)=BDATA(J) 
GO TO 98 
DO 99 I=1,NA 
BDATA(1)=BDATA(I )+DELTA/X 
CONVOLVE LOW PASS SMOOTHING WEIGHTS WITH ROUGH SPECTRUM 
NB=JX-N(K) 
GO TO 235 
IF (JCODE.EQ.1.AND.N(2).GT.0) GO TO 311 
IF(JCODE.LT.0) GO TO 311 
X=ALOG1O ( DELK/XORG)*3.0125 
TPP=BDATA(2)/YORG 
IF(TPP.LT.1.E-20) TPP=YORG 
Y=-0.5+AL0G10(TPP )*1. 35714 
CALL PLOT (XS Y, 3) 
D0 42 J=3,JX-1 
XJ=J-1 
X=AL0G10((XJ*DELK ) /XORG)*3.0125 
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TPP=BDATA(J)/YORG 
IF(TPP.LT.1.E-20) TPP=YORG 
Y=-0.5+AL0G10(TPP)*1. 35714 
42 CALL PLOT(X,Y,2) 
311 IF (KOUNT-2)78, 79,95 
78 KOUNT =2 
COMPUTE WEIGHTS FOR PREWHITENING FILTER AND STORE IN BDATA 
IF(N(1).LE.0) GO TO 95 
K=1 
98 NA=N(K)+1 
BDATA(1)=2.0* (CUT (K)+H(K) ) 
CENTER WEIGHT STORED IN LOCATION 1 
SUM=0.0 
PI=3. 1415926535898 
COMPUTE REMAINING WEIGHTS 
DO 40 I=2,NA 
P=I-1 
Q=1.0-16.*H(K)**2*P**2 
IF( ABS(Q).GT.0.0001) GO TO 11 
BDATA(I)= SIN(2.*PI*P* (CUT (K)+H(K)))/(4.0*P) 
GO TO 40 
11 BDATA(I)=(( COS(2.*PI*P*H(K)))/Q)*(( SIN(2.*PI*P*(CUT(K)+H(K)))) 
1 /(PI*P)) 
40 SUM=SUM+BDATA (I) 
CORRECT WEIGHTS FOR UNITY GAIN 
DELTA=1.-(BDATA(1)+2.*SUM) 
X=2*N(K)+1 
COMPUTE FINAL CORRECTED WEIGHTS ,HIGH PASS OR LOW PASS 
IF(K.EQ.2) GO TO 83 
DO 41 I=2,NA 
41 BDATA(I)= (BDATA(I)+ DELTA/X )*(-1.0) 
BDATA(1)=1.0-(BDATA(1)+DELTA/X ) 
NB=NP-N(K) 
CONVOLVE HIGH PASS WEIGHTS WITH ORIGINAL DATA,STORE IN ADATA 
235 DO 44 I=NA,NB 
IA=I-NA+1 
SUM=0.0 
DO 43 J=1,NA 
Jl=I+J-1 
J2=1-J+1 
43 SUM=SUM+BDATA(J)*(DATA(J1) + DATA(J2)) 
44 ADATA(IA)=SUM-BDATA(1)*DATA(I) 
IF(K.£Q.2) GO TO 236 
STORE PREWHITENED DATA IN DATA AND FILL IN ZEROES 
-NJ=NB-N(K ) 
DO 51 J=1,NJ 
51 DATA(J)=ADATA( J) 
NJ=NJ+1 
DO 54 J=NJ, INP 
54 DATA(J)=0.0 . 
PUT WEIGHTS IN PROPER ORDER -N.O.N AND STORE IN CDATA 
JJ=N(K)+1 
DO 52 J=1,JJ 
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491 


4911 


J1=J+N(K ) 
J2=J0= (0-1) 
CDATA(J1)=BDATA(J) 
CDATA(J2)=CDATA(J1) 
COMPUTE AMPLITUDE SPECTRUM OF PREWITENED DATA AND PLOT 
GO TO 53 
CONTINUE | 
STORE HIGH PASS WEIGHTS IN DATA AND AMP.SPECT.OF PREWHITENED 
DATA IN CDATA 
JJ=2*N(K )+1 
DO 102 J=1,JJ 
DATA(J)=CDATA(J) 
JJ=JJ+1 
DO 104 J=JJ,INP 
DATA(J)=0.0 
DO 103 J=1,JxX 
CDATA(J)=BDATA(J) 
COMPUTE AMPLITUDE SPECTRUM OF PREWHITENING FILTER 
KOUNT=3 
GO TO 53 
CONT INUE 
OUTPUT SMOOTHED SPECTRUM AND PLOT 
FIRST=(N(2)+1)*DELR 
FIRSTK=(N(2)+1)*DELK 
WRITE (IOUT, 237) FIRSTK 
FORMAT (38H FINAL SMOOTHED AND CORRECTED SPECT. ,1l6HFIRST FREQUENC 
1Y=,E8. 3) 
CALL PARSVL (ADATA,NP,INP,AVE,ANOIS ,N(2) ,ANORM ) 
WRITE (IOUT, 491) AVE ,ANOIS 
FORMAT(' VALUE OF AMP SPECT OF NOISE=',E10.4,' RMS WHITE NOISE LE 
1VEL=',£10.4) 
WRITE (IOUT, 81)(ADATA(J), J=2,1A) 
CALL POWWGT(B0,B81,FIRSTK,DELK,ADATA(2) ,IA-1,0,CTOFF1,CTOFF2) 
BOSAVE =B0 
BISAVE=B1 
CALCULATE RMS FOR WHITE NOISE WITH INTERCEPT (BO) 
AVE =B0 
AVE =AVE /ANORM 
AVE =AVE**2 
ANOISE=(AVE*INP)/(NP*INP ) 
ANOISE =SQRT (ANOISE ) 
WRITE (IOUT, 492)AVE,ANOISE 
FORMAT(' AVE= ',F8.3,' RMS WHITE NOISE = ',F8.3) 
IF(JCODE.LT.0) GO TO 486 
ENCODE(13,601,ITH) Bl 
FORMAT('SLOPE =',F6.3) 
CALL SYMBOL (7.5,8.5,.20,ITH,0.,13) 
ENCODE (19,602, 1TH)BO 
FORMAT('INTERCEPT =',£8.3) 
CALL SYMBOL (7.5,8.0,.20,ITH,0.,19) 
PLOT REGRESSION LINE FROM POWWGT 
CALL NEWPEN( 3) 
IPLPSS=0 
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603 FIRST1=FIRSTK 
C WRITE (IOUT, '(8HFIRST1= ,£11.4)' )FIRST1 
X=ALOG10(FIRST1/XORG)*3.0125 
AMP 1=BO*F IRST1**B1 
C WRITE (IOUT, '(6HAMP1= ,£11.4)')AMP1 
TPP=AMP1/YORG 
IF(TPP.LT.1.E-20) TPP=ABS(TPP) 
Y=-0.5+AL0G10(TPP )*1. 35714 
CALL PLOT(X,Y,+3) 
C CALL SYMBOL { X,Y,.1,32,0.0,-1) 
FIRST 1=((1A)*DELK ) 
C WRITE (IOUT, '(8HFIRST1= ,£11.4)')FIRST1 
X=ALOG10(FIRST1/XORG)*3.0125 
AMP 1=BO*F IRST1**B1 
C WRITE (IOUT, '(6HAMP1= ,£11.4)' )AMP1 
TPP=AMP 1/YORG 
IF(TPP.LT.1.E-20) TPP=ABS (TPP) 
Y=-0.5+AL0G10(TPP )*1. 35714 
CALL PLOT(X,Y,+2) 
CALL SYMBOL ( X,Y,.1,32,0.0,-2) 
CALL NEWPEN(1) 
PLOT SPECTRAL ESTIMATES 
X=ALOG10 (FIRSTK/XORG)*3.0125 
TPP=ADATA( 2) /YORG 
IF (TPP.EQ.0.0)TPP=1.E-15 
IF(TPP.LT.1.E-20) TPP=ABS (TPP) 
Y=-0.5+AL0G10(TPP)*1. 35714 
CALL PLOT(X,Y,+3) 
C ALL SYMBOL ( X,Y,.05,32,0.0,-1) 
DO 238 J=3,IA-1 
XJ=J-1+N(2) 
X=ALOG10( (XJ*DELK ) /XORG)*3.0125 
TPP=ADATA(J)/YORG 
IF (TPP.EQ.0.0)TPP=1.E-15 
IF(TPP.LT.1.E-20) TPP=ABS (TPP) 
Y=-0.5+ALOG10(TPP)*1. 35714 
238 CALL PLOT(X,Y,+2) 
C 238 CALL SYMBOL ( X,Y,.05,32,0.0,-2) 
CALL PLOT(0.0,0.0,+3) 
IF (IPLPSS.EQ.1) GO TO 330 
IF(LPASS.EQ.1) GO TO 330 
IPLPSS=1 
BO=XINTER 
B1=SLOPEX 
INCLUDE NEXT LINE IF A PLOT OF SPATIAL DOMAIN ESTIMATE IS DESIRED 
GO TO 603 
PLOT AXIS 
330 IF(ANORM.EQ.1.0) GO TO 333 
CALL AXES(0.,-.5, ‘NORMALIZED AMPLITUDE (IN KM)',27,9.5,90.0 
wes 4) O05 Os 0s— 1) 
GO TO 332 
333 CALL AXES(0.,-.5, AMPLITUDE (IN KILOMETERS)',25,9.5,90.0 
*,1.3514,10.0,0.0,-1) 


AND 
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C 
C 
C 


332 CALL AXES(0.,-0.5, 'FREQUENCY (CYCLES /KILOMETER)',27,12.2,0.0 
POOL 25 LOO NOnOee— 9) 
CALL SYMBOL (4.0,-1.4,.32,1TG,0.0, 30) 
AYORG=KYORG 
AXORG=KXORG 
00 362 J=1,8 
AA=J-1 
Y=(1.35714*AA)-0. 37 
FPN=AYORG+AA 
362 CALL NUMBER(-.25,Y,.16,FPN,90.0,-1) 
00 363 J=1,5 
AA=J-1 
X=(3.0125*AA )+0. 20 
FPN=AXORG+AA 
363 CALL NUMBER(X,-.25,.16,FPN,0.0,-1) 
DO 364 J=1,4 
IB=KXORG+J-1 
00 364 JC=2,9 
AC=JC 
AA=AC*10.0**IB 
X=AL0G10(AA/XORG)*3.0125 
CALL PLOT(X,-.5,3) 
364 CALL PLOT(X, -.46,2) 
00 365 J=1,7 
IB=KYORG+J-1 
DO 365 JC=2,9 
AC=JC 
AA=AC*10.0**IB 
Y=-0.5+AL0G10(AA/YORG)*1. 35714 
CALL PLOT(0.0,Y,3) 
365 CALL PLOT(-.04,Y,2) 
CALL PLOT(0.0,-1.0,3) 
CALL PLOT(0.0,-1.0,3) 


+ 


INTERACTIVE ROUTINE TO MODIFY INPUT STRING 


WRITE (6, 3641) 
3641 FORMAT(T70,' ARE YOU SATISFIED WITH THE DATA SET BOUNDARIES?'/ 
*T70,' ENTER E TO END RUN') 
READ(5,3652) RESPON 
IF(RESPON.EQ.'E') GO TO 487 
IF (RESPON.EQ.'Y') GO TO 3650 
WRITE (6, 3642) 
3642 FORMAT(T70,' WOULD YOU LIKE TO ELIMINATE THE ENTIRE SEGMENT?') 
READ(5, 3652) RESPON 
IF (RESPON.EQ.'N') GO TO 3644 
CALL ERASE 
GO TO 3658 
3644 WRITE (6, 3643) 
3643 FORMAT(T70,' ENTER FIRST,SECOND LOCATIONS(IN KM) FOR TRUNCATION' ) 
READ(5,*)FIR,SEC 
NF IR=FIR/GRIDKM 
IF(NFIR.LT.O) NFIR=0 
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NSEC=SEC/GRIDKM 
IF (NSEC.GT.NP) NSEC=NP 

C RECREATE DATA ARRAY 
DO 3645 ICHG=1,NP 
LOCCHG=ICHG+NFIR 
DATA( ICHG )=SDATA(LOCCHG ) 
IF (LOCCHG.EQ.NSEC) GO TO 3646 

3645 CONTINUE 
C INTERPOLATE NEW LAT & LON 
3646 SDIFF=FLOAT (NSEC )/FLOAT (NP) 

FDIFF =FLOAT (NFIR)/FLOAT (NP ) 
FLAT =DECDEG( ISDLAT ,SELAT ) 
FLON=DECDEG( ISDLNG,SELNG) 
SLAT =DECDEG(IDLAT ,EELAT ) 
SLON=DECDEG(IDLNG,EMLNG) 
FNLAT=((SLAT-FLAT )*FDIFF )+FLAT 
SNLAT=((SLAT -FLAT )*SDIFF )+FLAT 
FNLON=((SLON-FLON )*FDIFF )+FLON 
SNLON=((SLON-FLON )*SDIFF )+FLON 

C RECONSTRUCT INTO DEGREES AND MINUTES 
CALL DEGMIN(FNLAT, ISDLAT ,SELAT ) 
CALL DEGMIN(FNLON, ISDLNG,SELNG) 
CALL DEGMIN(SNLAT, IDLAT ,EELAT ) 
CALL DEGMIN(SNLON, IOLNG,EMLNG) 
NP=NSEC-NFIR 
CALL ERASE 
GO TO 2451 


INTERACTIVE ROUTINE fu MODIFY REGRESSION FIT TO SPECTRUM 


AND 


3650 WRITE (6, 3651) 
3651 FORMAT(T70,' ARE YOU SATISFIED WITH THE REGRESSION FIT?') 
READ(5, 3652 )RESPON 
3652 FORMAT (A1) 
IF(RESPON.EQ.'Y') GO TO 3656 
WRITE (6, 3653) 
3653 FORMAT(T70,' ENTER LOG(LOWER FREQUENCY) ,LOG(UPPER FREQUENCY )' ) 
READ(5,*) CTOFF1,CTOFF2 
CALL ERASE 
LPASS=1 
GO TO 4911 
3656 CTOFF 1=CT1SV 
CTOFF 2=CT 2SV 
CALL ERASE 
LPASS=0 
WRITE (14, 399)NP, IPNUM, ISDLAT ,SELAT , ISOLNG,SELNG, IDLAT,EELAT, 
LIDLNG,EMLNG, BISAVE , BOSAVE 
399 FORMAT(15,18,4(1X,14,F6.2),F6.3,£9.3) 
C PLOT PHASE SPECTRUM 
3658 IF(KPHA.NE.1) GO TO 581 
CALL PLOT(15.0,0.0,-3) 
X=ALOG10 (DELK /XORG)*3.0125 
Y=ABS (PHASE (2) /30.0) 
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IF(IPHA.GT.0) Y=Y*30.0/APHA 
IBB=4 
IF(Y.LT.8.0.AND.PHASE(2).GE.0.0) GO TO 552 
IF(Y.LT.8.0.AND.PHASE(2).LT.0.0) IBB=2 
IF(Y.GE.8.0) Y=8.5 
CALL SYMBOL (X,Y,.05,188,0.,-1) 
GO TO 553 
552 CALL PLOT(X,Y, 3) 
553 DO 571 J=3,JxX 
XJ=J-1 
X=AL0G10((XJ*DELK ) /XORG)*3.0125 
Y=ABS (PHASE (J )/30.0) 
IF(IPHA.GT.0) Y=¥*30.0/APHA 
IBB=4 
IF(Y.LT.8.0.AND.PHASE(J).GE.0.0) GO 
IF(Y.LT.8.0.AND.PHASE(J).LT.0.0) IBB 
IF(Y.GE.8.0) Y=8.5 
CALL SYMBOL (X,Y,.05,1BB,0.,-2) 
GO TO 571 
554 CALL PLOT(X,Y,2) 
571 CONTINUE 
486 CONTINUE 
COMPUTE AND PLOT SMOOTHED PHASE SPECTRUM 
IF(N(2).EQ.0) GO TO 579 
DO 572 I=NA,NB 
IA=I -NA+1 
SUM=0.0 
DO 573 J=1,NA 
Jl=I+J-1 
J2=I-J+1 
573 SUM=SUM+BDATA( J )* (PHASE (J1)+PHASE (J2) ) 
572 SMPH(IA)=SUM-BDATA(1)*PHASE (1) 
IF(IPHA.GT.0) GO TO 544 
WRITE (IOUT ,582) FIRSTK 
‘582 FORMAT(' SMOOTHED PHASE SPECTRUM(DEG)-FIRST FREQ.=',E8.3) 
GO TO 545 
544 WRITE(IOUT,546) FIRSTK 
546 FORMAT(' SMOOTHED PHASE SPECTRUM(DATA INT. )-FIRST FREQ.=',F8.5) 
545 WRITE( IOUT, 81) (SMPH(J),J=2,1A) 
IF(JCODE.LT.0) GO TO 581 
CALL NEWPEN(1) 
X=ALOG10(FIRSTK /XORG)*3.0125 
Y=ABS (SMPH(2)/30.0) 
IF(IPHA.GT.0) Y=Y¥*30.0/APHA 
IBB=4 
IF(Y.LT.8.0.AND.SMPH(2).GE.0.0) GO TO 548 
IF(Y.LT.8.0.AND.SMPH(2).LT.0.0) IBB=2 
IF(Y.GE.8.0) Y=8.5 
CALL SYMBOL (X,Y,.05,1BB,0.,-1) 
GO TO 549 
548 CALL PLOT(X,Y,3) 
549 DO 574 J=3,IA 
XJ=J-1+N(2) 


TO 554 
=2 
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X=AL0G10((XJ*DELK ) /XORG)*3. 0125 
Y=ABS (SMPH(J)/30.0) 
IF (IPHA.GT.0)Y=Y*30.0/APHA 
IBB=4 
IF (Y.LT.8.0.AND.SMPH(J).GE.0.0 
IF(Y.LT.8.0.AND.SMPH(J).LT.0.0 
IF(Y.GE.8.0) Y=8.5— 
CALL SYMBOL (X,Y,.05,IBB,0. ,-2) 
GO TO 574 
550 CALL PLOT(X,Y,2) 
574 CONTINUE 
PLOT PHASE AXIS 
579 IF(JCODE.LT.0) GO TO 581 
CALL NEWPEN(1) 
IF (IPHA.GT.0) GO TO 547 
CALL AXES(0.,0.,14HABS PHASE (DEG),14,6.,90.,1.,0.,30.,-1) 
GO TO 558 
547 CALL AXES(0.,0.,20HABS PHASE(DATA INT.),20,8.,90.,1.,0.,APHA, -1) 
558 CALL AXES(0.,0.,1TG,+30,12.2,0.,3.0125,10.,0.,-1) 
CALL SYMBOL (2.0,-.2,.14, 33HTRIANGLE INDICATES NEGATIVE PHASE,0.0, 
1 33) 
581 CONTINUE 
IF(N(2).EQ.0) GO TO 331 
COMPUTE AMPLITUDE SPECTRUM OF SMOOTHING FILTER 


wo 


GO TO 550 
IBB=2 


DELX= 0.01 

D0 239 J=1,51 
AJ=J-1 
XJ=AJd*DELX 
SUM=0.0 

DO 240 I=2,NA 
Al=I-1 


240 SUM=SUM+2*BDATA(I)* COS(2.*PI*AI*XJ) 
DATA(J)= SUM+BDATA(1) 
239 CDATA(J)=XJ 
WRITE(IOUT,241) CUT(2),H(2),N(2) 
241 FORMAT(' AMP.SPECT.OF SMOOTHING FILTER CUTOFF=',F5.4, 
We Ta Bein Gs Ue i avs 3053) 
IF (ILIST.NE.1) WRITE (IOUT, 242) (CDATA(I) ,DATA(I),1=1,51) 
242 FORMAT (8(F5.2,F8.3)) 
331 WRITE (IOUT, 243) JSET 
243 FORMAT(' END OF DATA SET NO.',I3) 
END FILE 12 
END FILE 12 
IF (JSET.EQ.NSETS) GO TO 244 
JSET=JSET +1 
IF(JCODE.LT.0) GO TO 245 
CALL PLOT(15.0,0.0,-3) 
GO TO 245 
244 IF(JCODE.LT.0) GO TO 487 
CALL PLOT(0.0,0.0,999) 
487 END FILE 14 
END FILE 14 
STOP 
END 
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THIS SUBROUTINE IS A SUPPLEMENT TO FFT1D PROGRAM 
AND DRAWS A PROFILE OF BATHYMETRIC DATA 
COMPILER (DIAG=3) 

SUBROUTINE BATHXS (DATA,NPTS ,GRIDKM ) 
DIMENSION DATA(1) 
CONVERT DEPTHS IN KILOMETERS TO METERS 
WRITE (IOUT, 3)(DATA(I),1=1,10) 
3 FORMAT(' ',10(F10.5,1X)) 
DO 5 I=1,NPTS 
5 DATA(I)=DATA(I)*1000. 
CALL PLOT(0.,10.,-3) 
FIND MAXIMUM AND MINIMUM DEPTH 
DE PMIN=0.0 
DE PMAX=0.0 
DO 10 I=1,NPTS 
IF(DATA(I).LT.DEPMIN) DEPMIN=DATA(I) 
10 IF(DATA(I).GT.DEPMAX ) DEPMAX=DATA(TI) 
ROUND TO NEAREST 100 
DE PMAX=(AINT (DEPMAX/100. )*100. )+100. 
DEPMIN=(AINT (DEPMIN/100. )*100. )-100. 
IF (GRIDKM.GT.0.0001) GO TO 11 
DE PMAX=-1.0 
DEPMIN=1.0 
GO TO 12 
11 IF(GRIDKM.GE.0.01) GO TO 12 
DE PMAX=-20.0 
DEPMIN=20.0 
12 YSCALE=-(DEPMAX-DEPMIN)/4. 

SETUP AXES 
ae AXES(0.,0.,' DEPTH (IN METERS)',18,4.0,90.,1.00,DEPMAX, YSCALE 

ae) 

XDIST=FLOAT (NPTS )*GRIDKM 
XSCALE =12.805/XDIST 
XINC=1.0 
55 XTICK=XSCALE 
56 IF(XTICK.GT.1.) GO TO 60 
XINC=2.*X INC 
XTICK=2.*XTICK 
GO T0 56 
60 IF(GRIDKM.GT.0.0001) GO TO 70 
CALL AXES(0.0,0.0,'DISTANCE (IN 0.1 METERS)',24,12.805,0.0,XTICK 
HOW oWai) 
G0 TO 90 
70 IF(GRIDKM.GT.0.01) GO TO 80 
CALL AXES(0.0,0.0, ‘DISTANCE (IN 100 METERS)',24,12.805,0.0,XTICK 
3560 3a 3) 
G0 TO 90 
80 CALL AXES(0.0,0.0,'DISTANCE (IN KILOMETERS )',24,12.805,0.0,XTICK 
*,0.0,XINC,-1) 
90 DO 100 I=1,NPTS 
YVAL =- ((DEPMAX-DATA(I))/YSCALE ) 
XVAL =(FLOAT (1 )*GRIDKM )*XSCALE 
IF(I.EQ.1) CALL PLOT(0.0,YVAL, 3) 
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C 


100 CALL PLOT(XVAL,YVAL, 2) 
CALL PLOT(0.,-10.,-3) 
CONVERT BACK TO KILOMETERS 
00 105 I=1,NPTS 

105 DATA(1I)=DATA(I)/1000. 
RETURN 
END 
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SUBROUTINE PRVOUT(NP,DEPTH,ITG,XINTER,SLOPE ,GRID) 
DIMENSION DEPTH(5000) 
COMMON /BOUND/1T, IPNUM, ISDLAT ,SELAT, ISDLNG,SELNG, IDLAT ,EELAT, 
1IDLNG,EMLNG 
CHARACTER*4 ITG(8) 
READER ROUTINE FOR DEPTH OUTPUT FROM PROVINCE PICKER 
READ(13,5,END=500) (ITG(I),1=1,6),SLOPE,XINTER,GRID,NP 
5 FORMAT (6A4,F7.3,1X,£9.3,1X,F7.4,15) 
NP=NP+1 
ITG(7)=' : 
ITG(8)=' : 
WRITE(6,55,END=500) (ITG(I),1=1,6) SLOPE ,XINTER,GRID,NP 
55 FORMAT (1T70,6A4,F7.3,1X,E10.5,F7.4,15) 
READ(13,50)1, IPNUM, ISDLAT ,SELAT, ISDLNG,SELNG, IDLAT ,EELAT, 
LIDLNG,EMLNG,RMS 
50 FORMAT (15,18,4(1X,14,F6.2),F 10.4) 
READ(13,6) (DEPTH(I),I=1,NP) 
6 FORMAT (10F 8.6) 
IF(NP.GT.2048) NP=2048 
READ(13,7) CHECK 
7 FORMAT (F8.2) 
IF (CHECK .NE.9.999999) PRINT 8 
8 FORMAT(' NINES RECORD DOES NOT CHECK‘) 
GO TO 900 
500 DEPTH(1)=9.999999 
900 RETURN 
END 
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C ROUTINE TO COMPUTE THE AVE VALUE OF THE LAST 10 PERCENT OF THE 
C HARMONICS(AVE) OF AN AMPLITUDE SPECTRUM(DATA) AND USE THIS VALUE WITH 
C PARSEVALS FORMULA TO COMPUTE THE RMS WHITE NOISE LEVEL(ANOIS). 
NP=THE NO.OF ORIGINAL PTS INPUT TO THE FFT USED TO COMPUTE DATA ARRAY 


C** 
C** 
C** 
C*x* 
Cx* 


SUBROUTINE PARSVL(DATA,NP,INP,AVE ,ANOIS ,NWGT ,ANORM ) 


INP=THE NO.OF POINTS TO A POWER OF 2 USED FOR THE FFT 
DATA= UNNORMALIZED FFT AMPLITUDE SPECTRUM WITH UNITS OF 


10 
11 


INPUT UNITS/CYCLE/DATA INTERVAL AND IS OF LENGTH(INP/2)+1-2*NWGT 
NWGT=LENGTH OF SMOOTHING FILTER YOU USED ON YOUR FFT -NWGT..0..+NWGT 


DIMENSION DATA(1) 

IF (ANORM.EQ.1.) GO TO 11 
DO 10 I=1,NP 

DATA(1 )=DATA(I) /ANORM 
JX=INP/2+1-2*NWGT 

JT =JX- (JX /10) 
X=0. 2* INP 

JX=X-2.0 

T=0.13*INP 

JT=T-2.0 

DIF =JX-JT+1 

ANP=NP 

AVE=0.0 

AVEPS=0.0 

00 1 I=JT,JX 
AVE=AVE+DATA(TI ) 
AVEPS=AVEPS+DATA(I )**2 
AVE =AVE /DIF 
AVEPS=AVEPS /DIF 


C PARSEVALS FORMULA IS MEAN SQ NOISE=(INP*AVEPS )/(NP*INP) 


ANOTS=SQRT (AVEPS /ANP ) 
IF (ANORM.EQ.1.) RETURN 
00 12 [=1,NP 


12 DATA(I )=DATA(1)*ANORM 


RETURN 
END 
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SUBROUTINE TSHIFT(M,DATA,A, DELF ) 
C*****THIS SUBROUTINE CONTAINS 26 STATEMENTS. 
C*****TSHIFT APPLIES EFFECTS OF TIME SHIFT ON REAL AND IMAGINARY PARTS OF THE 
C***x**FOURIER TRANSFORM. 
Cx****M -= POWER OF 2 WHICH IS EQUAL TO THE NUMBER OF REAL OR IMAGINARY PARTS 
(Chaat OF THE FOURIER TRANSFORM. 
C*****DATA -= SERIES ... REAL PARTS ARE ODD INDEX, IMAGINARY PARTS ARE EVEN. 
C*****xA -= YAXIS SHIFT, IE. TO CORRECT THE FOURIER TRANSFORM OF A FUNCTION 
WHICH 
Cae HAS BEEN SHIFTED A DATA INTERVALS IN THE +X DIRECTION,SIGN OF A IS 
+, 
C****x*DELF -= NORMALIZED FREQUENCY INCREMENT. 

DIMENSION DATA(1) 

M2=2%2%2M 

PT =3. 1415926536 

SFT =2. O*PI*A*DELF 

00 1 I=1,M2,2 

J=I+1 

TRE=DATA(I) 

TIM=DATA(J) 

K=1/2 

ARG=SFT *K 

CN=COS (ARG) 

SN=SIN (ARG) 

DATA(T )=TRE*CN-TIM*SN 

1 DATA(J)=TIM*CN+TRE*SN 
RETURN 
END 
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Appendix E 


In Chapter 6, a simple model of an anisotropic surface was created, 
and the variation of its spectral characteristics as a function of azi- 
muth was derived analytically and confirmed by empirical results. The 
surface was constructed by generating a random signal of known spectral 
properties and extending each value of the series into the second dimen- 


sion. In the resulting relationship, 


Ne a(8)°sP(9) 


The proportionality factor a(@) was shown to vary as a sinusoid and b(9) 
remained constant over all 0. Results from multibeam-sonar-derived 
bathymetry indicated that this simple model is adequate to describe some 
actual surfaces on ‘the earth, such as the Gorda Rise spreading center. 

Beyond this very simple model, a more elaborate surface can be gen- 
erated by summing several signals of differing characteristics at a 
variety of orientations. If one assumes that the simple corrugated 
surface presented in the elementary model of Chapter 6 is the result of 
a dominant, unidirectional process, then these composite surfaces would 
represent areas of the sea floor where several geological processes have 
been active, perhaps acting in different directions. The Gorda Rise, 
particularly at the ridge crest, is dominated by the processes associ- 
ated with crustal formation, and its relief does appear to conform to 
the one-process surface model. 

The relief of the Mendocino Fracture Zone shown in Mgures 4-7-10, 


although only represented by spectra from two azimuths, indicates a 
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clear difference in spectral slope (b(®)), which would not be expected 
in a simple surface dominated by one process. In fact, the main east- 
west trend of the fracture zone is produced by tectonic processes which 
are reflected in profiles collected in a north-south orientation. The 
various mass-wasting processes acting down this slope produce a differ- 
ent style of relief which is evident in profiles collected east-west. 
The bathymetry shown in Migure 6-10 from the continental slope might 
also be the results of two processes at work. This is reflected in the 
possible variation of b(®) as shown in Mgure 6-11. 

No attempt will be made to derive a general mathematical (geometri- 
cal) form for the azimuthally-dependent spectra of these multiple com- 
ponent surfaces. Such a treatment would be beyond the scope of this 
eoide. However, an insight into the nature of these surfaces can be 
gained by examining a few examples. All of the examples following are 
generated using the same algorithms. First, two random signals of spec- 
ified spectral characteristics (a and b) are generated using an inverse 
FFT method. One signal becomes the initial row, and the other the ini- 
tial column, of a 128 x 128 element matrix. Each element of the matrix 
is then computed by summing the appropriate row and column element. 
This method results in a surface which when sampled at O° azimuth, pro- 
duces one of the input signals displaced by a random constant. For the 
90° azimuth, the other input signal is sampled, similarly displaced. To 
study the spectra of the combined signals in other azimuths, the matri- 
ces were sampled and analyzed in an identical fashion to the bathymetry 
grids of Chapter 6. 

Four examples of two spectral component surfaces are shown in Fig- 


ures E-1 through E-4. A variety of combinations were used. Figure E-l 
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SLOPE OF SPECTRUM (b ) 


80 100 120 140 160 180 
AZIMUTH (degrees) 


20 40. 60 


Artificially generated surface composed of two identical 
orthogonal trends, both with spectral slope b= -1.5, and 
spectral intercept a= 1.0. Viewpoint is, from, the 
southwest. Below, the variation of parameters a and b are 
shown versus azimuth. 
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illustrates the characteristics of a surface in which two signals of 
identical spectral slope (b = -1.5) and intercept (a = 1.0) were com- 
bined in orthogonal directions. Arbitrarily assuming a viewpoint from 
the southwest for all examples, the input signals represent pure signals - 
in the east-west and north-south orientations. The resulting surface 
shows a clear northwest-southeast trend. Another realization might 
yield a northeast-southwest trend. With only a knowledge of the result- 
ing surface, an investigator might infer a single process acting at 45° 
or 135° azimuth. The trend actually results from the vector sum of two 
orthogonal processes. 

The distribution of spectral parameters with azimuth clearly 
reflects the departure of this surface from a one process model. The 
spectral slope parameter (b) does show the designated value of -1.5 at 
azimuths of 0°, 90° and 180°, as it must. The corresponding tatwecene™ 
(a) parameters aiso correspond to the input value of 1.0, indicating 
that these profiles represent uncontaminated samples of the input sig- 
nals. At intermediate azimuths, however, both parameters are consis- 
tently higher than those of the input signal. Whereas in the simple 
model of Chapter 6 the slope (b) parameter remains constant with azi- 
muth, the same parameter shows two clear maxima in the range 0° to 180°. 
The intercept (a) Parameter also shows two maxima, rather than the 
single maximum of the corrugated surface model. Figure E-2 presents a 
similar surface in which the intercept (a) parameter in the east-west 
direction has been increased to 1.2. The variation of spectral param- 
eters with azimuth is also very similar to that shown in Figure E-l, 


however the variation of both a and b is amplified. 
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SLOPE OF SPECTRUM (b) 
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AZIMUTH (degrees) 


Artificially generated surface composed of two orthogonal 


trends, both with spectral slopes of b= -1.5, and a=1.0 in 
the N-S direction, a=l1.2 in the E-W direction. Viewpoint 
is from the southwest. The variation of 4 and 6 parame- 
ters with azimuth are plotted below. : 
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In Figure E-3, the intercept terms (a) are again set equal, but the 
apectrall slopes (b) of the component signals differ. The north-south 
component has spectral slope b = -1.5, while the east-west component has 
a spectral slope b = -1.0. The artificial surface was designed to mimic 
the bathymetry of the Mendocino Fracture Zone. The north-south signal 
dominates the surface, as reflected in an east-west trending contour 
chart. The wavelength associated with the intercept term a, was arbi- 
trarily selected to correspond to a wavelength of two data points. 
Because of the higher angle of negative slope in the spectrum of the 
north-south profile, this trend contains higher amplitudes in all fre- 
quencies lower than one-cycle-per-two-data intervals. Only at the very 
southern limit of the surface, where the north-south signal is rela- 
tively constant, can the orthogonal trend be detected. 

At higher frequencies, the spectra of the two baupoueat signals 
cross over, and the east-west signal contains higher amplitudes and dou- 
inates the surface topography. This somewhat complicated set of circum- 
stances is not outside geological experience. For a terrestrial exam- 
ple, envision a long ridge of several kilometers width and perhaps one 
thousand meters height, oriented east-west. On the side of this ridge 
are a series of north-south trending streams and gullies with relief of 
tens-of-meters which shed the runoff from the ridge. If one were to 
travel in the north-south direction, the effect of the streams would be 
minimal, and the long wavelength shape of the ridge would present the 
only obstacle. If one were to travel in the east-west direction on the 
face of the ridge, the obstacles in the terrain would be dominated by 
the lower amplitude, but higher frequency streambeds. This difficulty 


of travel represents in a very direct sense the concept of surface 
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Figure E-3 Artificially generated surface composed of two orthogonal 
trends with identical intercepts of a= 1.0, but spectral 
slopes of b= -1.5 in the N-S direction and b= -1.0 in the 
E-W direction. Azimuthally dependent spectral parameters 
are plotted below. 
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roughness. The artificial surface shown in Figure E-3, as well as the 
Mendocino Fracture Zone shown in Chapter 4, are analogous to this hypo- 
thetical example. 

The computed spectral parameters shown below in Figure E-3 also 
reflect the complexities of this surface. The a(0) parameters appear to 
vary regularly with azimuth, although not following the cosine curve 
derived by simple regression, the variation occurs in spite of the fact 
that the input signals have identical a's. Only at exactly 90° azimuth, 
does the a parameter jump to the input value. The b(8) values also fol- 
low a complex pattern which is not described by the illustrated cosine 
curve. Extensive analytical geometry would be necessary to reach an 
understanding of these variations. 

A final artificial surface is presented as Migure E-4. In this 
case, both the spectral slopes (b) and intercepts (a) of the input sig- 
nals are d‘fferent The construction is identical to that shown in Fig- 
ure E-3, with the exception that the intercept (a) parameter in the 
east-west direction has been increased to 2.0. Between the input values 
at 0°, 90°, and 180°, the variation of the spectral parameters appears 
even more complicated than the results from Figure E-3. The combination 
of signals at oblique azimuths, or the combination of more than two sig- 
nals, would result in an even more complex pattern. 

With the insight gained by examining these artificially generated 
surfaces, a more complete analysis of the anisotropy of the Mendocino 
Fracture Zone can now be conducted. Figure E-5 presents a contour chart 
of the surface used in this analysis, which represents a subarea of the 
base chart shown in Figure 4-7. The digital data, collected by the SASS 


multibeam sonar system, is gridded at a spacing of 0.05 minutes of lati- 
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tude and longitude. Figure E-6 shows a graphic representation of this 
surface and its spectral parameters as a function of azimuth, in the 
same format as the artificially generated surfaces shown in Figures E-1 
through E-4. 

Compare the bathymetric surface shown in Figure E-6 to the artific- 
tially generated surface shown in Figures E-3 and E-4. The overall mor- 
phologies are quite similar, with a longer wavelength, higher amplitude 
component in the north-south direction relative to the orthogonal trend. 
This similarity in morphologies is also expressed as a similarity in the 
azimuthally dependent roughness models, although the parameters gener- 
ated from the bathymetric surface are somewhat noisier than the artifi- 
cially generated examples. In all cases, the slope parameter b is not 
constant with azimuth, as was the case for the surfaces studied in Chap- 
ter 6. Like the artificial surfaces of Figures E-3 and E-4, the spec- 
tral slope is approximately b = -1.5 in the 0° or 180° azimuth and 
approaches b = -1.0 for the 90° azimuth. The intercept parameter in 
Figure E-6 reaches its maximum at ® =90°, similar to the example in 
Figure E-4. In the case of the Mendocino Fracture Zone, this parameter 
is approximately doubled in the east-west direction over the north-south 
direction. This indicates that for wavelengths near 1 km, the Fracture 
Zone surface is twice as rough for ® = 90° as for 9 = 0°. In longer 
wavelengths, this relationship is reversed as evidenced by the much 
higher total relief in the north-south direction. Such a reversal is 
only possible for near orthogonal trends with different spectral slopes. 
Sinusoidal regression lines, which comprise the basic model of Chapter 
6, are included in Figure E-6 to emphasize how poorly this simple model 


describes the two-trend case. 
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Figure E-4 Artificially generated surface composed of orthogonal 
trends with spectral parameters b= -1.5, a= 1.0 in the N-S 
direction, and b= -1.0, a= 2.0 in the E-W direction. Azi- 
muthally dependent spectral parameters are plotted below. 
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Figure E-5 Contour chart of a segment of the Mendocino Fracture Zone. 
Automatically generated contours are based on SASS multibeam 


sonar data gridded at 0.05 minutes of latitude and 
longitude. 
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tude and longitude. Figure E-6 shows a graphic representation of this 
surface and its spectral parameters as a function of azimuth, in the 
same format as the artificially generated surfaces shown in Figures E-1 
through E-4. 

Compare the bathymetric surface shown in Figure E-6 to the artific- 
ially generated surface shown in Figures E-3 and E-4. The overall mor- 
phologies are quite similar, with a longer wavelength, higher amplitude 
component in the north-south direction relative to the orthogonal trend. 
This similarity in morphologies is also expressed as a similarity in the 
azimuthally dependent roughness models, although the parameters gener- 
ated from the bathymetric surface are somewhat noisier than the artific- 
jally generated examples. In all cases, the slope parameter b is not 
constant with azimuth, as was the case for the surfaces studied in Chap- 
ter 6. Like the artificial. surfaces of Figures E-3 and E-4, the spec- 
tral slope is approximately b = -1.5 in the 0° or 180° azimuth and 
approaches b = -1.0 for the 90° azimuth. The intercept parameter in 
Figure E-6 reaches its maximum at @ =90°, similar to the example in 
Figure E-4. In the case of the Mendocino Fracture Zone, this parameter 
is approximately doubled in the east-west direction over the north-south 
direction. This indicates that for wavelengths near 1 km, the Fracture 
Zone surface is twice as rough for 6 = 90° as for 8 = 0°. In longer 
wavelengths, this relationship is reversed as evidenced by the much 
higher total relief in the north-south direction. Such a reversal is 
only possible for near orthogonal trends with different spectral slopes. 
Sinusoidal regression lines, which comprise the basic model of Chapter 
6, are included in Figure E-6 to emphasize how poorly this simple model 


describes the two trend case. 
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Figure E-6 Graphic representation of the Mendocino Fracture Zone 
bathymetry shown in Figure E-5 is shown above (viewpoint 
from the southwest). Shown below are the spectral estimates 
versus azimuth for this surface. Compare this example to 
the artificially generated surfaces shown in Figures E-3 and 
E-4. 
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DESCRIPTION OF TERMS 


This informal glossary is included to aid readers who represent 
diverse backgrounds in geology, geophysics, acoustics, statistics, and 
other fields. The descriptions of terms, phrases, and acronyms included 
are intended to reflect their use within this report, rather than a 
general or rigorous definition. 


Amplitude 
The departure of a periodic function from its mean. 


Anisotropy 
Condition of having different properties in different directions. 


Band-1imited 
Condition in which a signal contains a finite range of frequency. 
With reference to Fourier analysis, this band ranges from the fun- 
damental frequency to the Nyquist frequency. 


Boxcar function 
A function which equals unity over some finite length and zero 
everywhere else. The multiplication of this function with an infi- 
nite series represents mathematically the sampling of a finite 
series from an infinite process. Sometimes called a rectangle 
function. 


Chi-square distribution 
The probability distribution for the estimation orror of the 
power spectrum. Due to the form of the distribution, the errors 
associated with the amplitude spectrum appear constant when plotted 
on log-transformed axes. 


Coherent signal 
A signal in which the phase relationship of the various frequency 
components is retained. 


Convolution — 
A mathematical operation which is equivalent to multiplication in 
the opposite transform domain. 


Deterministic model 
A numerical model in which the mathematical parameters represent 
specific measurable quantities of the process under study. 


Ensemble average 
The mean value of a group of repetitive samplings for a given sta- 
tistical measure. 


FFT 
Fast Fourier Transform. A numerical algorithm for estimating the 
Fourier transform with a greatly reduced number of operations. 
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Fractal dimension 
A topological term which refers to a dimension which which may be 
either integer or fractional. The fractal dimension (D) is related 
to the spectral slope (b) by b=-(5/2-D). Refer to Mandelbrot 
(1982) for a complete discussion. 


Functional form 
Refers to the type of function or functions used as a model for the 
distribution of data in a regression analysis. The term does not 
apply to the calculated parameters used to describe a specific data 
set. 


Fundamental frequency 
The lowest frequency treated in a Fourier analysis, corresponding 
to the inverse of the length of data. 


HEBBLE 
High Energy Benthic Boundary Layer Experiment 


Isotropy 
Condition of having the same properties in all directions. 


Leakage 
In spectral analysis, the transfer of energy from one frequency 
into other frequency bands. 


Linear-linear space 
A two-dimensional coordinate system in which both orthogonal axes 
represent simple evenly-spaced scales. Usually called a rec- 
tangular Cartesian coordinate system. 


Log-log space 
A two-dimensional coordinate system in which both orthogonal axes 
are scaled by evenly spacing the logarithm of the linear scale. In 
this study, all such transformations use base-ten logarithms. 


Magnetic anomaly ; 
The départure of the measured magnetic field from some low fre- 
quency model. 


Markov process 
A stochastic process in which the conditional probability state is 
unaffected by the historical state of the system. 


Multibeam sonar 
A bathymetric sounding system in which several discrete soundings 
of the sea floor can be derived by a single discharge of acoustic 
energy. 


NAVOCEANO 
United States Naval Oceanographic Office 


Non-parametric statistics 


Statistical theory in which the probability distribution of the 
underlying data is not assumed. 


216 


Nyquist frequency i 
The highest frequency treated in a Fourier analysis, corresponding 
to the inverse of twice the data spacing. 


Orthogonal 
An orientation in which all axes intersect perpendicularly. 


Phase 
The location of a periodic function along an axis relative to 
some arbitrary origin. 


Planetary Rossby waves 
A large-scale, stable wave motion in the global ocean. 


Power 
The squared amplitude of a periodic function. 


Power law 
A mathematical relationship of the form y=ax), This function 
plots linearly in log-log space. 


Prewhitening 
A technique which reduces the effect of spectral leakage in ana- 
lyzing non-white-noise signals. 


Process (geological) 
A natural continuing activity or function. 


Process (statistical) 
Any quantity which is defined in terms of its relationship to some 
independent variable, usually time or space. 


Provincing : 
Techniques which divide a sample space into discrete regions based 
on some predefined statistical property. 


Random walk 
A stochastic process of independent increments. 


Regression model 
A functional description of a relationship between a dependent 
variable and independent variable(s), usually derived by the method 
of least-squares. 


Round-off error 
The level of uncertainty in a data set due to the finite number of 
digits retained. 


SASS 
Sonar Array Sub-System. A multibeam sonar system operated by the 
U.S. Navy. 

SEABEAM 


A commercially available multibeam sonar system. 
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SEAMARC-1 
A commercially available deep-towed geophysical instrument package. 


Sedimentary process 
Relief-forming process which involves the transportation of 
suspended material. 


Sinc function 
The function y = sin(x)/x. 


Sinusoid 
A function having the form of a sine or cosine function, but with a 
variable phase value. 


Spatial frequency 
The inverse of wavelength. 


Spectral intercept 
The intersection of an amplitude spectrum with the amplitude axis. 
This term is specific to this report. 


Spectral slope 
The slope of an amplitude spectrum which has been plotted on log- 
log axes. This value becomes the exponent of spatial frequency 
following transformation to linear-linear space. This term is spe- 
cific to this report. 


Spreading center 
A region of the sea floor where recent crustal material is being 
formed. 


Stationarity 
Although rigorously defined in statistical theory, used in this 
report to imply relative constancy of a particular statistical 
parameter as a function of position; statistically homogeneous. 


Stereo-pair bottom photography 
A photogrammetric technique which allows the measurement of micro- 
relief on the sea floor by analyzing photographic images of a sur- 
face from offset viewpoints. 


Stochastic model 
A numerical model in which the mathematical parameters describe the 
process under study in terms of its random variability. 


Strike 
The bearing of the long axis of a linear trend. 


Tectonic process 
Relief-forming process which involves the deformation of the 
earth's crust. 


White noise 


A random series whose amplitude spectrum is constant with fre- 
quency. 
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